setup

# devtools::install_github("cran/drc", force = TRUE)
# devtools::install_github("cran/alr3", force = TRUE)
# devtools::install_github("cran/NRAIA", force = TRUE)
# devtools::install_github("cran/nlrwr", force = TRUE)
library(nlrwr)
Loading required package: alr3
Loading required package: car
Loading required package: carData
Loading required package: drc
Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:alr3’:

    forbes


'drc' has been loaded.

Please cite R and 'drc' if used for a publication,
for references type 'citation()' and 'citation('drc')'.


Attaching package: ‘drc’

The following objects are masked from ‘package:stats’:

    gaussian, getInitial

Loading required package: HydroMe
Loading required package: lattice
Loading required package: lmtest
Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading required package: NISTnls
Loading required package: nlme
Loading required package: nls2
Loading required package: proto
Loading required package: nlstools

'nlstools' has been loaded.

IMPORTANT NOTICE: Most nonlinear regression models and data set examples
related to predictive microbiolgy have been moved to the package 'nlsMicrobio'

Loading required package: NRAIA
Registered S3 method overwritten by 'NRAIA':
  method           from 
  plot.profile.nls stats

Attaching package: ‘NRAIA’

The following object is masked from ‘package:nlstools’:

    plotfit

Loading required package: sandwich

'nlrwr' has been loaded.

Please cite R and 'nlrwr' if used for a publication,
for references type 'citation()' and 'citation('nlrwr')'.

Attaching package: ‘nlrwr’

The following object is masked from ‘package:nlstools’:

    confint2

1 Introduction

1.1 A stock-recruitment model

## Fig. 1.1. 
plot(num.fish ~ spawn.biomass, 
     data = M.merluccius, 
     xlab = "Spawning biomass (1000 tonnes)",
     ylab = "Recruitment (million fish)")

1.2 Competition between plant biotypes

## Fig. 1.2.
plot(biomass ~ x, 
     data = RScompetition, 
     log = "", 
     xlab = Density ~ (plants/m^2), 
     ylab = Biomass ~ of ~ sensitive ~ biotype ~ (g/plant), 
     pch = as.numeric(as.factor(RScompetition$z)))

1.3 Grouped dose-response data

## Fig. 1.3.
xyplot(DryMatter ~ Dose | Herbicide,
       data = S.alba, 
       scales = list(x = list(log = TRUE)),
       ylab = "Dry matter (g/pot)",
       xlab = "Dose (g/ha)")

2 Getting Started

2.2 Getting started with nls()

2.2.1 Introducing the data example

L.minor
## Fig. 2.1.
plot(rate ~ conc, data = L.minor,
     ylab = "Uptake rate (weight/h)", 
     xlab = Substrate ~ concentration ~ (mmol ~ m^-3))

2.2.2 Model fitting

L.minor.m1 <- nls(rate ~ Vm * conc/(K + conc), 
                  data = L.minor, 
                  start = list(K = 20, Vm = 120), 
                  trace = TRUE)
624.3282    (1.32e+00): par = (20 120)
244.5460    (2.24e-01): par = (15.92382 124.5715)
234.5198    (2.91e-02): par = (17.25299 126.4388)
234.3595    (5.72e-03): par = (17.04442 125.9618)
234.3533    (1.11e-03): par = (17.08574 126.0467)
234.3531    (2.15e-04): par = (17.07774 126.0302)
234.3531    (4.19e-05): par = (17.0793 126.0334)
234.3531    (8.14e-06): par = (17.07899 126.0328)
deviance(L.minor.m1)
[1] 234.3531
d2 <- sum((L.minor$rate -predict(L.minor.m1))^2)
all.equal(deviance(L.minor.m1), d2)
[1] TRUE
logLik(L.minor.m1)
'log Lik.' -24.86106 (df=3)
coef(L.minor.m1)
        K        Vm 
 17.07899 126.03276 
summary(L.minor.m1)

Formula: rate ~ Vm * conc/(K + conc)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
K    17.079      2.953   5.784  0.00117 ** 
Vm  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 8.141e-06

2.2.3 Prediction

fitted(L.minor.m1)
[1]  18.06066  28.56474  38.52679  71.09461  78.03806  87.78424  91.62683 116.28685
attr(,"label")
[1] "Fitted values"
concVal <- with(L.minor, seq(min(conc), max(conc), length.out = 10))
predict(L.minor.m1, newdata = data.frame(conc = concVal))
 [1]  18.06066  75.09908  92.70509 101.26607 106.32776 109.67157 112.04518
 [8] 113.81733 115.19094 116.28685

2.2.4 Making plots

## Fig. 2.2.
plot(rate ~ conc, 
     data = L.minor,
     ylim = c(10, 130), 
     ylab = "Uptake rate (weight/h)", 
     xlab = Substrate ~ concentration ~ (mmol ~ m^-3)) 
lines(L.minor$conc, fitted(L.minor.m1))

## Fig. 2.3.
plot(rate ~ conc, 
     data = L.minor,
     ylim = c(10, 130), 
     ylab = "Uptake rate (weight/h)", 
     xlab = Substrate ~ concentration ~ (mmol ~ m^-3))

concVal <- with(L.minor, seq(min(conc), max(conc), length.out = 100))
lines(concVal, predict(L.minor.m1, newdata = data.frame(conc = concVal)))
abline(h = coef(L.minor.m1)[2], lty = 2)

2.2.5 Illustrating the estimation

L.minor.m1con <- nlsContourRSS(L.minor.m1)
0%100%
 RSS contour surface array returned 
## Fig. 2.4. 
par(pty = "s")
plot(L.minor.m1con, col = FALSE, nlev = 10)

2.3 Generalised linear models

## this code should be runned in the clean R session
## library(nlrwr) cause error 
# data(L.minor, package = "nlrwr")
# L.minor.m4 <- glm(rate ~ I(1/conc), 
#                   data = L.minor, 
#                   family = gaussian("inverse"))
# summary(L.minor.m4)

Excercises

## 2.1
y0_start = mean(RGRcurve[RGRcurve$Day == 0, "RGR"])
y1 = mean(RGRcurve[RGRcurve$Day == 1, "RGR"])
b_start = 1 / log(y1/y0_start)


e2.1 <- nls(RGR ~ y0 * exp(Day/b),
            data = RGRcurve,
            start = list(y0 = y0_start, b = b_start))

plot(RGRcurve$Day, RGRcurve$RGR, xlab = "Day", ylab = "RGR")
nd <- with(RGRcurve, seq(min(Day), max(Day), length.out = 101))
lines(nd, predict(e2.1, newdata = data.frame(Day = nd)), col = 2)

## 2.3
plotfit(L.minor.m1)

3 Starting Values and Self-starters

3.1 Finding starting values

3.1.1 Graphical exploration

## Fig. 3.1.
# install.packages("NISTnls")
data(Chwirut2, package = "NISTnls")
plot(y ~ x, 
     data = Chwirut2, 
     xlab = "Metal distance", 
     ylab = "Ultrasonic response")

## Fig. 3.2.
expFct <- function(x, beta1, beta2, beta3) {
  exp(-beta1 * x)/(beta2 + beta3 * x)
}

plot(y ~ x, 
     data = Chwirut2, 
     xlab = "Metal distance", 
     ylab = "Ultrasonic response",
     ylim = c(0, 100)
     )
curve(expFct(x, beta1 = 1, beta2 = 0.01, beta3 = 1), add = TRUE, col = 2)

## Fig. 3.3.
plot(y ~ x, data = Chwirut2, xlab = "Metal distance",
     ylab = "Ultrasonic response", ylim = c(0, 100))
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 1), add = TRUE, lty = 2,
      col = 2)
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 0.1), add = TRUE, lty = 3,
      col = 3)
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 0.01), add = TRUE, lty = 4,
      col = 4)
curve(expFct(x, beta1 = 0.2, beta2 = 0.01, beta3 = 0.01), add = TRUE, lty = 1,
      col = 6)

Chwirut2.m1 <- nls(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = list(beta1 = 0.2, beta2 = 0.01, beta3 = 0.01)
)
summary(Chwirut2.m1)

Formula: y ~ expFct(x, beta1, beta2, beta3)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 0.1665753  0.0383032   4.349 6.56e-05 ***
beta2 0.0051653  0.0006662   7.753 3.54e-10 ***
beta3 0.0121501  0.0015304   7.939 1.81e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.172 on 51 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 4.684e-06

3.1.2 Searching a grid

## bad starting example
# Chwirut2.m0 <- nls(
#   y ~ expFct(x, beta1, beta2, beta3),
#   data = Chwirut2,
#   start = list(beta1 = 1, beta2 = 0.01, beta3 = 1))
# summary(Chwirut2.m0)

## nlsLM may solve the bad starting
Chwirut2.m0 <- minpack.lm::nlsLM(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = list(beta1 = 1, beta2 = 0.01, beta3 = 1))
summary(Chwirut2.m0)

Formula: y ~ expFct(x, beta1, beta2, beta3)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 0.1665765  0.0383033   4.349 6.56e-05 ***
beta2 0.0051653  0.0006662   7.753 3.54e-10 ***
beta3 0.0121500  0.0015304   7.939 1.81e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.172 on 51 degrees of freedom

Number of iterations to convergence: 22 
Achieved convergence tolerance: 1.49e-08
grid.Chwirut2 <- expand.grid(beta1 = seq(0.1, 1, by = 0.1),
                             beta2 =  0.01,
                             beta3 = seq(0.1, 1, by = 0.1))
grid.Chwirut2
Chwirut2.m2a <- nls2(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = grid.Chwirut2,
  algorithm = "brute-force"
)
Chwirut2.m2a
Nonlinear regression model
  model: y ~ expFct(x, beta1, beta2, beta3)
   data: Chwirut2
beta1 beta2 beta3 
 0.10  0.01  0.10 
 residual sum-of-squares: 60696

Number of iterations to convergence: 100 
Achieved convergence tolerance: NA
nls(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = as.list(coef(Chwirut2.m2a)))  # starting value from Chwirut2.m2a
Nonlinear regression model
  model: y ~ expFct(x, beta1, beta2, beta3)
   data: Chwirut2
   beta1    beta2    beta3 
0.166578 0.005165 0.012150 
 residual sum-of-squares: 513

Number of iterations to convergence: 6 
Achieved convergence tolerance: 4.429e-06

3.2 Using self-starter functions

3.2.1 Built-in self-starter functions for nls()

L.minor.m2 <- nls(rate ~ SSmicmen(conc, Vm, K), data = L.minor)
summary(L.minor.m2)

Formula: rate ~ SSmicmen(conc, Vm, K)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
Vm  126.033      7.173  17.570 2.18e-06 ***
K    17.079      2.953   5.784  0.00117 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 2.478e-06

3.2.2 Defining a self-starter function for nls()

expModel <- function(predictor, b, y0) {
  y0 * exp(predictor/b)
}

expModelInit <- function(mCall, LHS, data) {
  xy <- sortedXyData(mCall[["predictor"]], LHS, data)
  lmFit <- lm(log(xy[, "y"]) ~ xy[, "x"])
  coefs <- coef(lmFit)
  y0 <- exp(coefs[1])
  b <- 1/coefs[2]
  value <- c(b, y0)
  names(value) <- mCall[c("b", "y0")]
  value
}

SSexp2 <- selfStart(expModel, expModelInit, c("b", "y0"))
with(RGRcurve, SSexp2(Day, 4, 0.2))
 [1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2568051
 [8] 0.2568051 0.2568051 0.2568051 0.2568051 0.3297443 0.3297443 0.3297443
[15] 0.3297443 0.3297443 0.3297443 0.4234000 0.4234000 0.4234000 0.4234000
[22] 0.4234000 0.4234000 0.5436564 0.5436564 0.5436564 0.5436564 0.5436564
[29] 0.5436564 0.8963378 0.8963378 0.8963378 0.8963378 0.8963378 0.8963378
[36] 1.4778112 1.4778112 1.4778112 1.4778112 1.4778112 1.4778112

Warning: selfStart caused an ERROR

## ERROR
# getInitial(RGR ~ SSexp2(Day, b, y0), data = RGRcurve)
#RGRcurve.m1 <- nls(RGR ~ SSexp2(Day, b, y0), data = RGRcurve)
#coef(RGRcurve.m1)

More examples on defining self-starter functions for nls() are found in Watkins and Venables (2006) or Venables and Ripley (2002a, pp. 216–217).

4 More on nls()

4.2 Supplying gradient information

4.2.1 Manual supply

MMfct1 <- function(conc, K, Vm) {
  numer <- Vm * conc
  denom <- K + conc
  mean <- numer/denom
  partialK <- -numer/(denom^2)
  partialVm <- mean/Vm
  attr(mean, "gradient") <- cbind(partialK, partialVm)
  return(mean)
}
MMfct1(1, 2, 3)
[1] 1
attr(,"gradient")
       partialK partialVm
[1,] -0.3333333 0.3333333
L.minor.mgr1 <- nls(rate ~ MMfct1(conc, K, Vm), 
                    data = L.minor, 
                    start = list(K = 1, Vm = 1),
                    trace = TRUE)
42634.85    (7.03e+00): par = (1 1)
31698.57    (9.47e+00): par = (389.6406 100.9883)
21357.28    (1.23e+01): par = (31.6085 46.75042)
18634.08    (1.03e+01): par = (27.26998 51.07311)
13989.24    (8.02e+00): par = (22.1864 59.6233)
7637.536    (5.64e+00): par = (18.35113 75.51722)
2003.136    (2.75e+00): par = (17.01927 100.4895)
234.3615    (6.31e-03): par = (17.1086 126.0567)
234.3533    (9.56e-04): par = (17.07326 126.021)
234.3531    (1.86e-04): par = (17.08017 126.0352)
234.3531    (3.61e-05): par = (17.07882 126.0324)
234.3531    (7.02e-06): par = (17.07909 126.0329)
summary(L.minor.mgr1)

Formula: rate ~ MMfct1(conc, K, Vm)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
K    17.079      2.953   5.784  0.00117 ** 
Vm  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 11 
Achieved convergence tolerance: 7.025e-06
## without derivative
nls(rate ~ Vm * conc/(K + conc), 
    data = L.minor, start = list(K = 1, Vm = 1), 
    trace = TRUE)
42634.85    (7.03e+00): par = (1 1)
31698.57    (9.47e+00): par = (389.6406 100.9883)
21357.27    (1.23e+01): par = (31.60848 46.75041)
18634.07    (1.03e+01): par = (27.26996 51.07311)
13989.24    (8.02e+00): par = (22.18639 59.6233)
7637.536    (5.64e+00): par = (18.35113 75.51722)
2003.136    (2.75e+00): par = (17.01927 100.4895)
234.3615    (6.31e-03): par = (17.1086 126.0567)
234.3533    (9.56e-04): par = (17.07326 126.021)
234.3531    (1.86e-04): par = (17.08017 126.0352)
234.3531    (3.61e-05): par = (17.07882 126.0324)
234.3531    (7.03e-06): par = (17.07909 126.0329)
Nonlinear regression model
  model: rate ~ Vm * conc/(K + conc)
   data: L.minor
     K     Vm 
 17.08 126.03 
 residual sum-of-squares: 234.4

Number of iterations to convergence: 11 
Achieved convergence tolerance: 7.029e-06

4.2.2 Automatic supply

MMfct2 <- deriv(~ Vm * conc / (K + conc), 
                c("K", "Vm"), 
                function(conc, K, Vm) {})
MMfct2
function (conc, K, Vm) 
{
    .expr1 <- Vm * conc
    .expr2 <- K + conc
    .value <- .expr1/.expr2
    .grad <- array(0, c(length(.value), 2L), list(NULL, c("K", 
        "Vm")))
    .grad[, "K"] <- -(.expr1/.expr2^2)
    .grad[, "Vm"] <- conc/.expr2
    attr(.value, "gradient") <- .grad
    .value
}
MMfct2(1, 2, 3)
[1] 1
attr(,"gradient")
              K        Vm
[1,] -0.3333333 0.3333333
L.minor.mgr2 <- nls(
  rate ~ MMfct2(conc, + K, Vm), 
  data = L.minor, 
  start = list(K = 20, Vm = 120))
summary(L.minor.mgr2)

Formula: rate ~ MMfct2(conc, +K, Vm)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
K    17.079      2.953   5.784  0.00117 ** 
Vm  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 8.148e-06

Ref: Venables와 Ripley(2002a, pp. 214–216)

4.3 Conditionally linear parameters

4.3.1 nls() using the “plinear” algorithm

L.minor.m3 <- nls(rate ~ conc/(K + conc), 
                  data = L.minor, 
                  algorithm = "plinear",
                  start = list(K = 20))
summary(L.minor.m3)

Formula: rate ~ conc/(K + conc)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
K      17.079      2.953   5.784  0.00117 ** 
.lin  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 2.143e-06

4.3.2 A pedestrian approach

## Fig. 4.1.
data(segreg)
plot(C ~ Temp, data = segreg, 
     xlab = Mean ~ temperature ~ (degree ~ F),
     ylab = "Energy consumption (kWh)")

profRSS1 <- function(gamma) {
  deviance(lm(C ~ pmax(0, Temp - gamma), data = segreg))
}
profRSS2 <- Vectorize(profRSS1, "gamma")
l <- lm( C ~ pmax(0, Temp - 40), data = segreg)
summary(l)

Call:
lm(formula = C ~ pmax(0, Temp - 40), data = segreg)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.472  -2.155  -0.148   2.678  11.775 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        74.43931    1.16428  63.936  < 2e-16 ***
pmax(0, Temp - 40)  0.53635    0.06226   8.614 2.27e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.312 on 37 degrees of freedom
Multiple R-squared:  0.6673,    Adjusted R-squared:  0.6583 
F-statistic: 74.21 on 1 and 37 DF,  p-value: 2.268e-10
plot(segreg$Temp, segreg$C)
lines(segreg$Temp, predict(l), col = 2)

## Fig. 4.2. 
plot(profRSS2(Temp) ~ Temp, 
     data = segreg,
     type = "l", 
     xlab = expression(gamma),
     ylab = "Profile RSS")

## 좀더 부드럽게 만듬
plot(profRSS2(Temp) ~ Temp, 
     data = data.frame(Temp = seq(min(segreg$Temp), 
                                  max(segreg$Temp), 
                                  length.out = 1001)),
     type = "l", 
     xlab = expression(gamma),
     ylab = "Profile RSS")

4.4 Fitting models with several predictor variables

4.4.1 Two-dimensional predictor

RScompetition.m1 <- nls(
  biomass ~ a/(1 + b * (x + c * z)),
  data = RScompetition,
  start = list(a = 20, b = 1, c = 1)
)
summary(RScompetition.m1)

Formula: biomass ~ a/(1 + b * (x + c * z))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a  25.9144     2.3393  11.078 1.42e-14 ***
b   0.1776     0.0355   5.004 8.68e-06 ***
c   1.1349     0.2964   3.829 0.000387 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.673 on 46 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 1.437e-06
virDensity <- with(RScompetition, x + coef(RScompetition.m1)[3] * z)
virDenVal <- seq(0, max(virDensity), length.out = 100)
biomassVal <- predict(RScompetition.m1, data.frame(x = virDenVal, z = 0))
## Fig. 4.3. 
plot(biomassVal ~ virDenVal, type = "l", 
     ylab = "Biomass of sensitive biotype (g/plant)",
     xlab = Virtual ~ density ~ (plants/m^2))
with(RScompetition, points(biomass ~ virDensity))

4.4.2 General least-squares minimisation

## Fig. 4.4. 
par(pty = "s")
plot(Q ~ I, data = IQsig)
theta <- 0:360 * (pi/180)
lines(cos(theta), sin(theta))

IQsig.m1 <- nls(
  ~ ((I - I0)^2 - 2 * gamma * sin(phi) * (I - I0) * (Q - Q0) + 
       gamma * gamma * (Q - Q0)^2) - (rho * gamma * cos(phi))^2, 
  data = IQsig, 
  start = list(I0 = -0.005, gamma = 1, phi = -0.005, Q0 = -0.005, rho = 1))
summary(IQsig.m1)

Formula: 0 ~ ((I - I0)^2 - 2 * gamma * sin(phi) * (I - I0) * (Q - Q0) + 
    gamma * gamma * (Q - Q0)^2) - (rho * gamma * cos(phi))^2

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
I0    -0.002259   0.002058  -1.097    0.278    
gamma  1.007053   0.003719 270.764  < 2e-16 ***
phi   -0.053540   0.004977 -10.758 5.02e-14 ***
Q0    -0.002481   0.002278  -1.089    0.282    
rho    0.974295   0.002540 383.553  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01917 on 45 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 4.699e-07

5 Model Diagnostics

5.2 Checking the mean structure

5.2.1 Plot of the fitted regression curve

## Fig. 5.1.
plot(p ~ T, data = vapCO, log = "y",
     xlab = "Temperature (K)",
     ylab = "Pressure (Pa)")

vapCO.m1 <- nls(
  log(p) ~ A - B/(C + T), 
  data = vapCO,
  start = list(A = 10, B = 100, C = -10))
summary(vapCO.m1)

Formula: log(p) ~ A - B/(C + T)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
A  19.0923     0.2312   82.56  < 2e-16 ***
B 497.8246    27.7002   17.97 4.84e-10 ***
C -16.1516     1.5437  -10.46 2.19e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09951 on 12 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 3.585e-06
## Fig. 5.2.
plot(p ~ T, data = vapCO, log = "y",
     xlab = "Temperature (K)",
     ylab = "Pressure (Pa)")
lines(vapCO$T, exp(fitted(vapCO.m1)), col = 2)

## Fig. 5.3. 
plot(weight ~ conc, data = lettuce,
     xlab = "Concentration (mg/l)",
     ylab = "Biomass (g)", log = "x")
Warning in xy.coords(x, y, xlabel, ylabel, log) :
  2 x values <= 0 omitted from logarithmic plot

5.2.2 Residual plots

## Fig. 5.4.
plot(fitted(vapCO.m1), residuals(vapCO.m1),
     xlab = "Fitted Values", ylab = "Residuals")
abline(a = 0, b = 0)

## Fig. 5.5. 
plot(rootl ~ conc, data = ryegrass,
     xlab = "Concentration (mM)",
     ylab = "Root length (cm)")

ryegrass.m1 <- lm(rootl ~ as.factor(conc), data = ryegrass)
summary(ryegrass.m1)

Call:
lm(formula = rootl ~ as.factor(conc), data = ryegrass)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.08968 -0.17121  0.01577  0.16847  1.21865 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          7.74935    0.22535  34.388  < 2e-16 ***
as.factor(conc)0.94 -0.07606    0.39032  -0.195  0.84780    
as.factor(conc)1.88 -1.33479    0.39032  -3.420  0.00327 ** 
as.factor(conc)3.75 -4.73466    0.39032 -12.130 8.53e-10 ***
as.factor(conc)7.5  -6.71542    0.39032 -17.205 3.45e-12 ***
as.factor(conc)15   -7.07018    0.39032 -18.114 1.50e-12 ***
as.factor(conc)30   -7.44601    0.39032 -19.077 6.47e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.552 on 17 degrees of freedom
Multiple R-squared:  0.9791,    Adjusted R-squared:  0.9718 
F-statistic:   133 on 6 and 17 DF,  p-value: 2.488e-13
ryegrass.m2 <- nls(
  rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e)))),
  data = ryegrass,
  start = list(b = 1, c = 0.6, d = 8, e = 3)
)
summary(ryegrass.m2)

Formula: rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e))))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
b   2.9822     0.4584   6.505 2.43e-06 ***
c   0.4814     0.2095   2.298   0.0325 *  
d   7.7930     0.1897  41.073  < 2e-16 ***
e   3.0580     0.1858  16.456 4.31e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5196 on 20 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 6.33e-07
anova(ryegrass.m2, ryegrass.m1)
Analysis of Variance Table

Model 1: rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e))))
Model 2: rootl ~ as.factor(conc)
  Res.Df Res.Sum Sq Df  Sum Sq F value Pr(>F)
1     20     5.4002                          
2     17     5.1799  3 0.22035  0.2411 0.8665
Q <- -2 * (logLik(ryegrass.m2) - logLik(ryegrass.m1))
df.Q <- df.residual(ryegrass.m2) - df.residual(ryegrass.m1)
1 - pchisq(Q, df.Q)
'log Lik.' 0.8012878 (df=5)

5.3 Variance homogeneity

5.3.2 Levene’s test

with(ryegrass, leveneTest(rootl, as.factor(conc)))
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  6  1.9266 0.1344
      17               
## Fig. 5.6. 
plot(fitted(vapCO.m1), abs(residuals(vapCO.m1)),
     xlab = "Fitted values", ylab = "Absolute residuals")

5.4 Normal distribution

## Fig. 5.7. 
plot(vapCO.m1)

5.4.1 QQ plot

## Fig. 5.8. 
standardRes <- residuals(ryegrass.m2)/summary(ryegrass.m2)$sigma

par(pty = "s")
qqnorm(standardRes, main = "")
abline(a = 0, b = 1)

5.4.2 Shapiro-Wilk test

shapiro.test(standardRes)

    Shapiro-Wilk normality test

data:  standardRes
W = 0.98232, p-value = 0.9345

5.5 Independence

## Fig. 5.9. 
plot(residuals(vapCO.m1), c(residuals(vapCO.m1)[-1], NA),
     xlab = "Residuals", ylab = "Lagged residuals")

6 Remedies for Model Violations

6.1 Variance modelling

6.1.1 Power-of-the-mean variance model

## Fig. 6.1. 
plot(RGR ~ Day, data = RGRcurve,
     xlab = "Time (days)", ylab = "Relative growth rate (%)")

y0_start = mean(RGRcurve[RGRcurve$Day == 0, "RGR"])
y1 = mean(RGRcurve[RGRcurve$Day == 1, "RGR"])
b_start = 1 / log(y1/y0_start)

RGRcurve.m2 <- gnls(
  RGR ~  y0 * exp(Day/b), 
  data = RGRcurve,
  start = list(y0 = y0_start, b = b_start),
  weights = varPower()
)
summary(RGRcurve.m2)
Generalized nonlinear least squares fit
  Model: RGR ~ y0 * exp(Day/b) 
  Data: RGRcurve 

Variance function:
 Structure: Power of variance covariate
 Formula: ~fitted(.) 
 Parameter estimates:
   power 
1.012821 

Coefficients:

 Correlation: 
  y0   
b 0.797

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-1.8322566 -0.6356287 -0.1612768  0.6476686  1.8176689 

Residual standard error: 0.1672069 
Degrees of freedom: 41 total; 39 residual
RGRcurve.m2g <- gnls(
  RGR ~  y0 * exp(Day/b), 
  data = RGRcurve,
  start = list(y0 = y0_start, b = b_start)
  # weights = varPower() # without weights
)
summary(RGRcurve.m2g )
Generalized nonlinear least squares fit
  Model: RGR ~ y0 * exp(Day/b) 
  Data: RGRcurve 

Coefficients:

 Correlation: 
  y0   
b 0.962

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-2.4946822 -0.5595546 -0.1065972  0.3415913  2.3870757 

Residual standard error: 0.1009825 
Degrees of freedom: 41 total; 39 residual
AIC(RGRcurve.m2, RGRcurve.m2g )

6.2 Transformations

6.2.1 Transform-both-sides approach

## Fig. 6.2. 
plot(recruits ~ spawners, 
     data = sockeye[-12,],
     xlab = "Number of spawners (thousands)",
     ylab = "Number of recruits (thousands)")

일단 nls를 이용하여 fitting을 한 후

sockeye.m1 <- nls(
  recruits ~ beta1 * spawners * exp(-beta2 * spawners),
  data = sockeye[-12, ],
  start = list(beta1 = 2, beta2 = 0.001)
)
summary(sockeye.m1)

Formula: recruits ~ beta1 * spawners * exp(-beta2 * spawners)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 3.9495589  1.0154151   3.890 0.000658 ***
beta2 0.0008526  0.0003570   2.388 0.024812 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 566.1 on 25 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 3.796e-07
## Fig. 6.3. 
plot(fitted(sockeye.m1), abs(residuals(sockeye.m1)),
     xlab = "Fitted values", ylab = "Absolute residuals")

6.2.2 Finding an appropriate transformation

boxcox.nls를 이용하여 양쪽을 변화시키며 피팅을 한다.

## Fig. 6.4.
sockeye.m2 <- boxcox.nls(sockeye.m1)

bcSummary(sockeye.m2)

Estimated lambda: -0.2 
Confidence interval for lambda: [-0.89, 0.52] 
coef(summary(sockeye.m1))
          Estimate   Std. Error  t value    Pr(>|t|)
beta1 3.9495588725 1.0154150955 3.889600 0.000657508
beta2 0.0008525523 0.0003570088 2.388043 0.024811780
coef(summary(sockeye.m2))
          Estimate   Std. Error  t value     Pr(>|t|)
beta1 3.7767882170 0.6941022768 5.441256 1.195243e-05
beta2 0.0009539343 0.0003061565 3.115839 4.563494e-03

6.3 Sandwich estimators

이 내용은 정확히 이해하지는 못했음.

vcov(sockeye.m1)
             beta1        beta2
beta1 1.0310678161 3.436055e-04
beta2 0.0003436055 1.274553e-07
sandwich(sockeye.m1)
             beta1        beta2
beta1 0.6082739015 2.458975e-04
beta2 0.0002458975 1.162302e-07
coeftest(sockeye.m1)

t test of coefficients:

        Estimate Std. Error t value  Pr(>|t|)    
beta1 3.94955887 1.01541510  3.8896 0.0006575 ***
beta2 0.00085255 0.00035701  2.3880 0.0248118 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
coeftest(sockeye.m1, vcov = sandwich)

t test of coefficients:

        Estimate Std. Error t value  Pr(>|t|)    
beta1 3.94955887 0.77991916  5.0641 3.158e-05 ***
beta2 0.00085255 0.00034093  2.5007   0.01931 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

6.4 Weighting

6.4.1 Decline in nitrogen content in soil

exp1
## Fig. 6.5. 
plot(Nremaining ~ time, data = exp1,
     xlab = "Time (years)", ylab = "Nitroge content (%)")

exp1.m1 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp1)
exp1.m2 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp1,
               weights = norep/(stdev * stdev))
weights(exp1.m2)
[1] 0.4403928 0.4403928 0.1096776 0.1430179 1.1291355 0.0867694 1.8896447 0.3667661
coef(summary(exp1.m1))
     Estimate Std. Error    t value    Pr(>|t|)
a1 49.6945735  6.1714388  8.0523481 0.001291449
a2  0.2182533  0.3162611  0.6901049 0.528080839
b1 32.2059455  6.1631160  5.2255946 0.006402929
b2 -3.1242900  0.4881789 -6.3998879 0.003061147
coef(summary(exp1.m2))
      Estimate Std. Error     t value    Pr(>|t|)
a1 48.95580529  5.8454395  8.37504275 0.001111764
a2 -0.02288415  0.2953566 -0.07747973 0.941962763
b1 30.67388498  6.4550371  4.75193009 0.008958175
b2 -3.15057479  0.5205589 -6.05229290 0.003760909
exp2
exp2.m1 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp2)
exp2.m2 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp2,
               weights = norep/(stdev^2))
## Fig. 6.6. 
plot(Nremaining ~ time, data = exp2,
     xlab = "Time (years)", ylab = "Nitrogen content (%)")
timeVal <- with(exp2, seq(min(time), max(time), length.out = 101))
lines(timeVal, predict(exp2.m1, newdata = data.frame(time = timeVal)), 
      lty = 2, col = 2)
lines(timeVal, predict(exp2.m2, newdata = data.frame(time = timeVal)), 
      lty = 3, col = 3)

coef(summary(exp2.m1))
     Estimate Std. Error     t value    Pr(>|t|)
a1 37.1828528  7.3946625   5.0283367 0.007342087
a2  0.3984263  0.4676840   0.8519136 0.442262753
b1 56.9940958  7.2995778   7.8078620 0.001451984
b2 -2.4593036  0.2429193 -10.1239548 0.000535818
coef(summary(exp2.m2))
     Estimate Std. Error     t value     Pr(>|t|)
a1 28.9158392  8.3745109   3.4528392 2.598700e-02
a2  0.3128281  0.8122243   0.3851499 7.197304e-01
b1 59.7777230  3.7595533  15.9002197 9.144765e-05
b2 -2.3567707  0.1311759 -17.9664898 5.641332e-05
summary(exp1.m1)

Formula: Nremaining ~ SSbiexp(time, a1, a2, b1, b2)

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
a1  49.6946     6.1714   8.052  0.00129 **
a2   0.2183     0.3163   0.690  0.52808   
b1  32.2059     6.1631   5.226  0.00640 **
b2  -3.1243     0.4882  -6.400  0.00306 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.844 on 4 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 5.364e-06
summary(exp1.m2)

Formula: Nremaining ~ SSbiexp(time, a1, a2, b1, b2)

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
a1 48.95581    5.84544   8.375  0.00111 **
a2 -0.02288    0.29536  -0.077  0.94196   
b1 30.67388    6.45504   4.752  0.00896 **
b2 -3.15057    0.52056  -6.052  0.00376 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.785 on 4 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 9.037e-06

7 Uncertainty, Hypothesis Testing, and Model Selection

7.1 Profile likelihood

L.minor.m1 <- update(L.minor.m1, trace = FALSE)
L.minor.m1pro <- profile(L.minor.m1)
## Fig. 7.1. 
plot(L.minor.m1pro)

## Fig. 7.2. 
plot(L.minor.m1pro, absVal = FALSE)

confint(L.minor.m1)
Waiting for profiling to be done...
        2.5%     97.5%
K   11.45142  24.86684
Vm 110.64232 143.75648

7.2 Bootstrap

set.seed(222)
L.minor.m1boot <- nlsBoot(L.minor.m1)
summary(L.minor.m1boot)

------
Bootstrap statistics
    Estimate Std. error
K   17.37534   2.730157
Vm 126.64954   6.522141

------
Median of bootstrap estimates and percentile confidence intervals
      Median      2.5%     97.5%
K   17.25307  12.30746  23.20204
Vm 126.08120 115.00875 139.46932
## Fig. 7.3. 
par(mfrow = c(1, 2))
qqnorm(L.minor.m1boot$coefboot[, 1], main = "K")
qqnorm(L.minor.m1boot$coefboot[, 2], main = "Vm")

7.3 Wald confidence intervals

confint2(L.minor.m1)
        2.5 %    97.5 %
K    9.853828  24.30416
Vm 108.480813 143.58470
confint2(L.minor.m1, level = 0.99)
       0.5 %    99.5 %
K   6.131815  28.02617
Vm 99.439005 152.62651

7.4 Estimating derived parameters

deltaMethod(L.minor.m1, "Vm/(4*K)")
           Estimate     SE  2.5 % 97.5 %
Vm/(4 * K)   1.8449 0.2363 1.3817  2.308

7.5 Nested models

secalonic
## Fig. 7.4.
plot(rootl ~ dose, data = secalonic,
     xlab = "Dose (mM)", ylab = "Root length (cm)")

7.5.1 Using t-tests

secalonic.m1 <- nls(rootl ~ SSfpl(dose, a, b, c, d), data = secalonic)
summary(secalonic.m1)

Formula: rootl ~ SSfpl(dose, a, b, c, d)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a 6.053612   0.395467  15.308 0.000606 ***
b 0.353944   0.194089   1.824 0.165722    
c 0.075188   0.005911  12.721 0.001048 ** 
d 0.029350   0.006621   4.433 0.021333 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2028 on 3 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 8.102e-06

7.5.2 Using F-tests

secalonic.m2 <- nls(rootl ~ SSlogis(dose, a, c, d), data = secalonic)
summary(secalonic.m2)

Formula: rootl ~ SSlogis(dose, a, c, d)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a  6.344853   0.564690  11.236 0.000357 ***
c  0.076959   0.009397   8.190 0.001211 ** 
d -0.035740   0.007244  -4.933 0.007854 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2546 on 4 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 4.999e-07
anova(secalonic.m2, secalonic.m1)
Analysis of Variance Table

Model 1: rootl ~ SSlogis(dose, a, c, d)
Model 2: rootl ~ SSfpl(dose, a, b, c, d)
  Res.Df Res.Sum Sq Df  Sum Sq F value Pr(>F)
1      4    0.25924                          
2      3    0.12341  1 0.13582  3.3016 0.1668

7.6 Non-nested models

M.merluccius
plot(num.fish ~ spawn.biomass, data = M.merluccius)

4가지 조금씩 다른 모델을 비교

M.merluccius.bh <- nls(
  num.fish ~ spawn.biomass * alpha/(1 + spawn.biomass/k),
  data = M.merluccius, start = list(alpha = 5, k = 50)
)

M.merluccius.de <- nls(
  num.fish ~ spawn.biomass * alpha * (1 - c * spawn.biomass/k)^(1/c),
  data = M.merluccius, start = list(alpha = 4.4, k = 106, c = 0.86)
)

M.merluccius.ri <- nls(
  num.fish ~ spawn.biomass * alpha * exp(-spawn.biomass/k),
  data = M.merluccius, start = list(alpha = 5, k = 50)
)

M.merluccius.sh <- nls(
  num.fish ~ spawn.biomass * alpha/(1 + (spawn.biomass/k)^c),
  data = M.merluccius, start = list(alpha = 3.87, k = 61.72, c = 2.25),
  control = nls.control(maxiter = 100)
)

Residual standard error를 계산하면,

summary(M.merluccius.bh)$sigma
[1] 14.69956
summary(M.merluccius.de)$sigma
[1] 14.78346
summary(M.merluccius.ri)$sigma
[1] 14.30425
summary(M.merluccius.sh)$sigma
[1] 14.56758

AIC를 비교해보면,

AIC(M.merluccius.bh)
[1] 127.0562
AIC(M.merluccius.de)
[1] 128.0263
AIC(M.merluccius.ri)
[1] 126.2383
AIC(M.merluccius.sh)
[1] 127.585

Burnham and Anderson (2002, pp. 70–72)은 한 모델이 다른 모델에 비하여 절대적으로 우위에 있기 위해서는 AIC의 차이가 10 이상의 차이가 필요하다는 경험적 법칙을 제시하였다. 따라서 4 모델 사이에서 큰 차이는 없다.

8 Grouped Data

8.1 Fitting grouped data models

Puromycin
## Fig. 8.1.
xyplot(rate ~ conc | state, data = Puromycin,
       xlab = "Substrate concentration (ppm)",
       ylab = "Reaction rates\n(counts/min/min)")

8.1.1 Using nls()

Puromycin.m1 <- nls(
  rate ~ Vm[state] * conc/(K[state] + conc),
  data = Puromycin,
  start = list(K = c(0.1, 0.1), Vm = c(200, 200))
)
summary(Puromycin.m1)

Formula: rate ~ Vm[state] * conc/(K[state] + conc)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
K1  6.412e-02  7.877e-03   8.141 1.29e-07 ***
K2  4.771e-02  8.281e-03   5.761 1.50e-05 ***
Vm1 2.127e+02  6.608e+00  32.185  < 2e-16 ***
Vm2 1.603e+02  6.896e+00  23.242 2.04e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.4 on 19 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 1.568e-06

8.1.2 Using gnls()

Puromycin.m2 <- gnls(
  rate ~ Vm * conc/(K + conc),
  data = Puromycin,
  start = list(Vm = c(200, 200), K = c(0.1, 0.1)),
  params = list(Vm ~ state - 1, K ~ state - 1)
)
summary(Puromycin.m2)
Generalized nonlinear least squares fit
  Model: rate ~ Vm * conc/(K + conc) 
  Data: Puromycin 

Coefficients:

 Correlation: 
                  Vm.sttt Vm.sttn K.sttt
Vm.stateuntreated 0.000                 
K.statetreated    0.765   0.000         
K.stateuntreated  0.000   0.777   0.000 

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-1.32632770 -0.52742178 -0.06890298  0.51295804  2.44557119 

Residual standard error: 10.40003 
Degrees of freedom: 23 total; 19 residual

8.1.3 Using nlsList()

Puromycin.m3 <- nlsList(
  rate ~ SSmicmen(conc, Vm, K) | state,
  data = Puromycin
)
summary(Puromycin.m3)
Call:
  Model: rate ~ SSmicmen(conc, Vm, K) | state 
   Data: Puromycin 

Coefficients:
   Vm 
          Estimate Std. Error  t value     Pr(>|t|)
treated   212.6837   6.608092 32.18534 3.241160e-11
untreated 160.2801   6.896019 23.24241 1.384618e-09
   K 
            Estimate  Std. Error  t value     Pr(>|t|)
treated   0.06412123 0.007876786 8.140532 1.565136e-05
untreated 0.04770829 0.008281171 5.761056 1.727047e-04

Residual standard error: 10.40003 on 19 degrees of freedom
Puromycin2 <- groupedData(rate ~ conc | state, data = Puromycin)
Puromycin.m4 <- nlsList(
  rate ~ SSmicmen(conc, a, b), data = Puromycin2
)
summary(Puromycin.m4)
Call:
  Model: rate ~ SSmicmen(conc, a, b) | state 
   Data: Puromycin2 

Coefficients:
   a 
          Estimate Std. Error  t value     Pr(>|t|)
untreated 160.2801   6.896019 23.24241 1.384618e-09
treated   212.6837   6.608092 32.18534 3.241160e-11
   b 
            Estimate  Std. Error  t value     Pr(>|t|)
untreated 0.04770829 0.008281171 5.761056 1.727047e-04
treated   0.06412123 0.007876786 8.140532 1.565136e-05

Residual standard error: 10.40003 on 19 degrees of freedom

8.2 Model reduction and parameter models

8.2.1 Comparison of entire groups

Puromycin.m5 <- nls(
  rate ~ Vm * conc/(K + conc), data = Puromycin,
  start = list(K = 0.1, Vm = 200))
summary(Puromycin.m5)

Formula: rate ~ Vm * conc/(K + conc)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K    0.06039    0.01077   5.608 1.45e-05 ***
Vm 190.80620    8.76458  21.770 6.84e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18.61 on 21 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 6.259e-06
anova(Puromycin.m5, Puromycin.m1)
Analysis of Variance Table

Model 1: rate ~ Vm * conc/(K + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
  Res.Df Res.Sum Sq Df Sum Sq F value    Pr(>F)    
1     21     7276.5                                
2     19     2055.1  2 5221.5  24.138 6.075e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
AIC(Puromycin.m5, Puromycin.m1)

8.2.2 Comparison of specific parameters

Puromycin.m6 <- nls(
  rate ~ Vm * conc/(K[state] + conc), 
  data = Puromycin, 
  start = list(K = c(0.1, 0.1), Vm = 200))
summary(Puromycin.m6)

Formula: rate ~ Vm * conc/(K[state] + conc)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K1   0.04801    0.00874   5.493 2.24e-05 ***
K2   0.08865    0.01448   6.122 5.54e-06 ***
Vm 194.10328    7.40464  26.214  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.52 on 20 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 8.957e-06
anova(Puromycin.m6, Puromycin.m1)
Analysis of Variance Table

Model 1: rate ~ Vm * conc/(K[state] + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
  Res.Df Res.Sum Sq Df Sum Sq F value    Pr(>F)    
1     20     4815.4                                
2     19     2055.1  1 2760.3  25.521 7.082e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Puromycin.m7 <- nls(
  rate ~ Vm[state] * conc/(K + conc), 
  data = Puromycin, start = list(K = 0.1, Vm = c(200, 200))
)
summary(Puromycin.m7)

Formula: rate ~ Vm[state] * conc/(K + conc)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
K     0.05797    0.00591   9.809 4.37e-09 ***
Vm1 208.63004    5.80399  35.946  < 2e-16 ***
Vm2 166.60408    5.80743  28.688  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.59 on 20 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 1.364e-06
anova(Puromycin.m7, Puromycin.m1)
Analysis of Variance Table

Model 1: rate ~ Vm[state] * conc/(K + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
  Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1     20     2240.9                         
2     19     2055.1  1 185.84  1.7182 0.2056

AIC로 비교해 보아도 Puromycin.m7이 더욱 좋다.

AIC(Puromycin.m1, Puromycin.m5, Puromycin.m6, Puromycin.m7)

8.3 Common control

G.aparine
## Fig. 8.2. 
xyplot(drymatter ~ dose | as.factor(treatment),
       data = G.aparine,
       xlab = "Dose (g/ha)", ylab = "Dry matter (mg/pot)")

이미 데이터에서 treatment 0가 없어졌음. 아마도 데이터를 업데이트 한 것으로 추측됨.

G.aparine$treatment2 <- factor(G.aparine$treatment)
#levels(G.aparine$treatment2) <- c("0", "1", "2")
G.aparine
G.aparine.m1 <- nls(
  drymatter ~ c[treatment2] + (d - c[treatment2])/(1 + exp(b[treatment2] * (log(dose) - log(e[treatment2])))),
  data = G.aparine, 
  start = list(b = c(2, 2), c = c(500, 100), d = 1000, e = c(50, 100)))
summary(G.aparine.m1)

Formula: drymatter ~ c[treatment2] + (d - c[treatment2])/(1 + exp(b[treatment2] * 
    (log(dose) - log(e[treatment2]))))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b1   1.6129     0.3462   4.659 5.33e-06 ***
b2   1.7510     0.2303   7.603 7.07e-13 ***
c1 509.5035    23.4379  21.738  < 2e-16 ***
c2 151.9164    27.2114   5.583 6.56e-08 ***
d  984.8883    12.7957  76.970  < 2e-16 ***
e1  50.8000     7.8465   6.474 5.58e-10 ***
e2  93.4466     8.1626  11.448  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 111.4 on 233 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 9.945e-06

8.4 Prediction

concValues <- with(Puromycin, seq(min(conc), max(conc, length.out = 11)))
concValues
 [1]  0.02  1.02  2.02  3.02  4.02  5.02  6.02  7.02  8.02  9.02 10.02
stateVal1 <- levels(Puromycin$state)
stateVal1
[1] "treated"   "untreated"
csValues1 <- expand.grid(conc = concValues, state = stateVal1)
csValues1
predict(Puromycin.m7, newdata = csValues1)
 [1]  53.51423 197.41021 202.80963 204.70062 205.66419 206.24825 206.64012
 [8] 206.92127 207.13280 207.29773 207.42993  42.73445 157.64435 161.95611
[15] 163.46619 164.23566 164.70207 165.01500 165.23951 165.40844 165.54015
[22] 165.64572
stateVal2 <- factor("untreated", levels = c("treated", "untreated"))
stateVal2
[1] untreated
Levels: treated untreated
csValues2 <- data.frame(conc = concValues, state = stateVal2)
csValues2
predict(Puromycin.m1, newdata = csValues2)
 [1]  47.3444 153.1183 156.5819 157.7874 158.4002 158.7711 159.0198 159.1981
 [9] 159.3322 159.4367 159.5205

8.5 Nonlinear mixed models

vinclozolin
## Fig. 8.3.
xyplot(effect ~ conc | exper, data = vinclozolin,
       xlab = expression(paste("Concentration (", mu, "M)")),
       ylab = "Luminescence (LU")

LL3.formula <- effect ~ d/(1 + exp(b * (log(conc) - log(e))))

vinclozolin.e1.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10509,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e2.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10821,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e3.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10828,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e4.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10904,  
  start = list(b = 1, d = 2700, e = 0.03))
vinclozolin.e5.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 11023,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e6.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 11106,  
  start = list(b = 0.5, d = 2600, e = 0.02))
## Fig. 8.4. 
#par(pty = "s")
plot(effect ~ conc, data = vinclozolin,
     pch = as.numeric(exper), 
     log = "x",
     xlim = c(1e-04, 10), 
     xlab = expression(paste("Concentration(", mu, "M)")), 
     ylab = "Luminescence(LU)")
concVec <- exp(seq(log(1e-04), log(10), length.out = 50))
lines(concVec, predict(vinclozolin.e1.m, data.frame(conc = concVec)), lty = 2)
lines(concVec, predict(vinclozolin.e2.m, data.frame(conc = concVec)), lty = 3)
lines(concVec, predict(vinclozolin.e3.m, data.frame(conc = concVec)), lty = 4)
lines(concVec, predict(vinclozolin.e4.m, data.frame(conc = concVec)), lty = 5)
lines(concVec, predict(vinclozolin.e5.m, data.frame(conc = concVec)), lty = 6)
lines(concVec, predict(vinclozolin.e6.m, data.frame(conc = concVec)), lty = "3313")

vinclozolin.m1 <- nlme(
  effect ~ d/(1 + exp(b * (log(conc) - log(e)))), 
  data = vinclozolin,
  fixed = list(b ~ 1, d ~ 1, e ~ 1), 
  random = d ~ 1 | exper, 
  start = c(1, 1000, 0.1))
summary(vinclozolin.m1)
Nonlinear mixed-effects model fit by maximum likelihood
  Model: effect ~ d/(1 + exp(b * (log(conc) - log(e)))) 
  Data: vinclozolin 

Random effects:
 Formula: d ~ 1 | exper
               d Residual
StdDev: 695.3775 150.9623

Fixed effects:  list(b ~ 1, d ~ 1, e ~ 1) 
 Correlation: 
  b      d     
d -0.051       
e  0.529 -0.145

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.83842036 -0.58871949 -0.06161384  0.54652838  2.07902663 

Number of Observations: 53
Number of Groups: 6 
## Fig. 8.5.
#par(pty = "s")
plot(effect ~ conc, data = vinclozolin,
     pch = as.numeric(exper), 
     log = "x",
     xlim = c(1e-04, 10), 
     xlab = expression(paste("Concentration(", mu, "M)")), 
     ylab = "Luminescence(LU)")
concVec <- exp(seq(log(1e-04), log(10), length.out = 50))
lines(concVec, predict(vinclozolin.e1.m, data.frame(conc = concVec)), lty = 2)
lines(concVec, predict(vinclozolin.e2.m, data.frame(conc = concVec)), lty = 3)
lines(concVec, predict(vinclozolin.e3.m, data.frame(conc = concVec)), lty = 4)
lines(concVec, predict(vinclozolin.e4.m, data.frame(conc = concVec)), lty = 5)
lines(concVec, predict(vinclozolin.e5.m, data.frame(conc = concVec)), lty = 6)
lines(concVec, predict(vinclozolin.e6.m, data.frame(conc = concVec)), lty = "3313")
lines(concVec, 
      predict(vinclozolin.m1, newdata = data.frame(conc = concVec), level = 0), 
      lty = 1, lwd = 3)

—-

LS0tCnRpdGxlOiAiTm9ubGluZWFyIFJlZ3Jlc3Npb24gd2l0aCBSIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgoKIyBzZXR1cAoKYGBge3J9CiMgZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJjcmFuL2RyYyIsIGZvcmNlID0gVFJVRSkKIyBkZXZ0b29sczo6aW5zdGFsbF9naXRodWIoImNyYW4vYWxyMyIsIGZvcmNlID0gVFJVRSkKIyBkZXZ0b29sczo6aW5zdGFsbF9naXRodWIoImNyYW4vTlJBSUEiLCBmb3JjZSA9IFRSVUUpCiMgZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJjcmFuL25scndyIiwgZm9yY2UgPSBUUlVFKQpsaWJyYXJ5KG5scndyKQpgYGAKCgojIDEgSW50cm9kdWN0aW9uCgojIyAxLjEgQSBzdG9jay1yZWNydWl0bWVudCBtb2RlbAoKYGBge3J9CiMjIEZpZy4gMS4xLiAKcGxvdChudW0uZmlzaCB+IHNwYXduLmJpb21hc3MsIAogICAgIGRhdGEgPSBNLm1lcmx1Y2NpdXMsIAogICAgIHhsYWIgPSAiU3Bhd25pbmcgYmlvbWFzcyAoMTAwMCB0b25uZXMpIiwKICAgICB5bGFiID0gIlJlY3J1aXRtZW50IChtaWxsaW9uIGZpc2gpIikKYGBgCgoKIyMgMS4yIENvbXBldGl0aW9uIGJldHdlZW4gcGxhbnQgYmlvdHlwZXMKCmBgYHtyfQojIyBGaWcuIDEuMi4KcGxvdChiaW9tYXNzIH4geCwgCiAgICAgZGF0YSA9IFJTY29tcGV0aXRpb24sIAogICAgIGxvZyA9ICIiLCAKICAgICB4bGFiID0gRGVuc2l0eSB+IChwbGFudHMvbV4yKSwgCiAgICAgeWxhYiA9IEJpb21hc3MgfiBvZiB+IHNlbnNpdGl2ZSB+IGJpb3R5cGUgfiAoZy9wbGFudCksIAogICAgIHBjaCA9IGFzLm51bWVyaWMoYXMuZmFjdG9yKFJTY29tcGV0aXRpb24keikpKQpgYGAKCgojIyAxLjMgR3JvdXBlZCBkb3NlLXJlc3BvbnNlIGRhdGEKCmBgYHtyfQojIyBGaWcuIDEuMy4KeHlwbG90KERyeU1hdHRlciB+IERvc2UgfCBIZXJiaWNpZGUsCiAgICAgICBkYXRhID0gUy5hbGJhLCAKICAgICAgIHNjYWxlcyA9IGxpc3QoeCA9IGxpc3QobG9nID0gVFJVRSkpLAogICAgICAgeWxhYiA9ICJEcnkgbWF0dGVyIChnL3BvdCkiLAogICAgICAgeGxhYiA9ICJEb3NlIChnL2hhKSIpCmBgYAoKCiMgMiBHZXR0aW5nIFN0YXJ0ZWQKCiMjIDIuMiBHZXR0aW5nIHN0YXJ0ZWQgd2l0aCBubHMoKQoKIyMjIDIuMi4xIEludHJvZHVjaW5nIHRoZSBkYXRhIGV4YW1wbGUKCmBgYHtyfQpMLm1pbm9yCmBgYAoKCmBgYHtyfQojIyBGaWcuIDIuMS4KcGxvdChyYXRlIH4gY29uYywgZGF0YSA9IEwubWlub3IsCiAgICAgeWxhYiA9ICJVcHRha2UgcmF0ZSAod2VpZ2h0L2gpIiwgCiAgICAgeGxhYiA9IFN1YnN0cmF0ZSB+IGNvbmNlbnRyYXRpb24gfiAobW1vbCB+IG1eLTMpKQpgYGAKCgojIyMgMi4yLjIgTW9kZWwgZml0dGluZwoKYGBge3J9CkwubWlub3IubTEgPC0gbmxzKHJhdGUgfiBWbSAqIGNvbmMvKEsgKyBjb25jKSwgCiAgICAgICAgICAgICAgICAgIGRhdGEgPSBMLm1pbm9yLCAKICAgICAgICAgICAgICAgICAgc3RhcnQgPSBsaXN0KEsgPSAyMCwgVm0gPSAxMjApLCAKICAgICAgICAgICAgICAgICAgdHJhY2UgPSBUUlVFKQpgYGAKCgpgYGB7cn0KZGV2aWFuY2UoTC5taW5vci5tMSkKYGBgCgoKYGBge3J9CmQyIDwtIHN1bSgoTC5taW5vciRyYXRlIC1wcmVkaWN0KEwubWlub3IubTEpKV4yKQphbGwuZXF1YWwoZGV2aWFuY2UoTC5taW5vci5tMSksIGQyKQpgYGAKCgpgYGB7cn0KbG9nTGlrKEwubWlub3IubTEpCmBgYAoKCmBgYHtyfQpjb2VmKEwubWlub3IubTEpCmBgYAoKCmBgYHtyfQpzdW1tYXJ5KEwubWlub3IubTEpCmBgYAoKCiMjIyAyLjIuMyBQcmVkaWN0aW9uCgpgYGB7cn0KZml0dGVkKEwubWlub3IubTEpCmBgYAoKCmBgYHtyfQpjb25jVmFsIDwtIHdpdGgoTC5taW5vciwgc2VxKG1pbihjb25jKSwgbWF4KGNvbmMpLCBsZW5ndGgub3V0ID0gMTApKQpwcmVkaWN0KEwubWlub3IubTEsIG5ld2RhdGEgPSBkYXRhLmZyYW1lKGNvbmMgPSBjb25jVmFsKSkKYGBgCgoKIyMjIDIuMi40IE1ha2luZyBwbG90cwoKYGBge3J9CiMjIEZpZy4gMi4yLgpwbG90KHJhdGUgfiBjb25jLCAKICAgICBkYXRhID0gTC5taW5vciwKICAgICB5bGltID0gYygxMCwgMTMwKSwgCiAgICAgeWxhYiA9ICJVcHRha2UgcmF0ZSAod2VpZ2h0L2gpIiwgCiAgICAgeGxhYiA9IFN1YnN0cmF0ZSB+IGNvbmNlbnRyYXRpb24gfiAobW1vbCB+IG1eLTMpKSAKbGluZXMoTC5taW5vciRjb25jLCBmaXR0ZWQoTC5taW5vci5tMSkpCmBgYAoKCmBgYHtyfQojIyBGaWcuIDIuMy4KcGxvdChyYXRlIH4gY29uYywgCiAgICAgZGF0YSA9IEwubWlub3IsCiAgICAgeWxpbSA9IGMoMTAsIDEzMCksIAogICAgIHlsYWIgPSAiVXB0YWtlIHJhdGUgKHdlaWdodC9oKSIsIAogICAgIHhsYWIgPSBTdWJzdHJhdGUgfiBjb25jZW50cmF0aW9uIH4gKG1tb2wgfiBtXi0zKSkKCmNvbmNWYWwgPC0gd2l0aChMLm1pbm9yLCBzZXEobWluKGNvbmMpLCBtYXgoY29uYyksIGxlbmd0aC5vdXQgPSAxMDApKQpsaW5lcyhjb25jVmFsLCBwcmVkaWN0KEwubWlub3IubTEsIG5ld2RhdGEgPSBkYXRhLmZyYW1lKGNvbmMgPSBjb25jVmFsKSkpCmFibGluZShoID0gY29lZihMLm1pbm9yLm0xKVsyXSwgbHR5ID0gMikKYGBgCgoKIyMjIDIuMi41IElsbHVzdHJhdGluZyB0aGUgZXN0aW1hdGlvbgoKYGBge3J9CkwubWlub3IubTFjb24gPC0gbmxzQ29udG91clJTUyhMLm1pbm9yLm0xKQpgYGAKCgpgYGB7cn0KIyMgRmlnLiAyLjQuIApwYXIocHR5ID0gInMiKQpwbG90KEwubWlub3IubTFjb24sIGNvbCA9IEZBTFNFLCBubGV2ID0gMTApCmBgYAoKCiMjIDIuMyBHZW5lcmFsaXNlZCBsaW5lYXIgbW9kZWxzCgpgYGB7cn0KIyMgdGhpcyBjb2RlIHNob3VsZCBiZSBydW5uZWQgaW4gdGhlIGNsZWFuIFIgc2Vzc2lvbgojIyBsaWJyYXJ5KG5scndyKSBjYXVzZSBlcnJvciAKIyBkYXRhKEwubWlub3IsIHBhY2thZ2UgPSAibmxyd3IiKQojIEwubWlub3IubTQgPC0gZ2xtKHJhdGUgfiBJKDEvY29uYyksIAojICAgICAgICAgICAgICAgICAgIGRhdGEgPSBMLm1pbm9yLCAKIyAgICAgICAgICAgICAgICAgICBmYW1pbHkgPSBnYXVzc2lhbigiaW52ZXJzZSIpKQojIHN1bW1hcnkoTC5taW5vci5tNCkKYGBgCgoKIyMgRXhjZXJjaXNlcwoKYGBge3J9CiMjIDIuMQp5MF9zdGFydCA9IG1lYW4oUkdSY3VydmVbUkdSY3VydmUkRGF5ID09IDAsICJSR1IiXSkKeTEgPSBtZWFuKFJHUmN1cnZlW1JHUmN1cnZlJERheSA9PSAxLCAiUkdSIl0pCmJfc3RhcnQgPSAxIC8gbG9nKHkxL3kwX3N0YXJ0KQoKCmUyLjEgPC0gbmxzKFJHUiB+IHkwICogZXhwKERheS9iKSwKICAgICAgICAgICAgZGF0YSA9IFJHUmN1cnZlLAogICAgICAgICAgICBzdGFydCA9IGxpc3QoeTAgPSB5MF9zdGFydCwgYiA9IGJfc3RhcnQpKQoKcGxvdChSR1JjdXJ2ZSREYXksIFJHUmN1cnZlJFJHUiwgeGxhYiA9ICJEYXkiLCB5bGFiID0gIlJHUiIpCm5kIDwtIHdpdGgoUkdSY3VydmUsIHNlcShtaW4oRGF5KSwgbWF4KERheSksIGxlbmd0aC5vdXQgPSAxMDEpKQpsaW5lcyhuZCwgcHJlZGljdChlMi4xLCBuZXdkYXRhID0gZGF0YS5mcmFtZShEYXkgPSBuZCkpLCBjb2wgPSAyKQpgYGAKCgpgYGB7cn0KIyMgMi4zCnBsb3RmaXQoTC5taW5vci5tMSkKYGBgCgoKIyAzIFN0YXJ0aW5nIFZhbHVlcyBhbmQgU2VsZi1zdGFydGVycwoKIyMgMy4xIEZpbmRpbmcgc3RhcnRpbmcgdmFsdWVzCgojIyMgMy4xLjEgR3JhcGhpY2FsIGV4cGxvcmF0aW9uCgpgYGB7cn0KIyMgRmlnLiAzLjEuCiMgaW5zdGFsbC5wYWNrYWdlcygiTklTVG5scyIpCmRhdGEoQ2h3aXJ1dDIsIHBhY2thZ2UgPSAiTklTVG5scyIpCnBsb3QoeSB+IHgsIAogICAgIGRhdGEgPSBDaHdpcnV0MiwgCiAgICAgeGxhYiA9ICJNZXRhbCBkaXN0YW5jZSIsIAogICAgIHlsYWIgPSAiVWx0cmFzb25pYyByZXNwb25zZSIpCmBgYAoKCmBgYHtyfQojIyBGaWcuIDMuMi4KZXhwRmN0IDwtIGZ1bmN0aW9uKHgsIGJldGExLCBiZXRhMiwgYmV0YTMpIHsKICBleHAoLWJldGExICogeCkvKGJldGEyICsgYmV0YTMgKiB4KQp9CgpwbG90KHkgfiB4LCAKICAgICBkYXRhID0gQ2h3aXJ1dDIsIAogICAgIHhsYWIgPSAiTWV0YWwgZGlzdGFuY2UiLCAKICAgICB5bGFiID0gIlVsdHJhc29uaWMgcmVzcG9uc2UiLAogICAgIHlsaW0gPSBjKDAsIDEwMCkKICAgICApCmN1cnZlKGV4cEZjdCh4LCBiZXRhMSA9IDEsIGJldGEyID0gMC4wMSwgYmV0YTMgPSAxKSwgYWRkID0gVFJVRSwgY29sID0gMikKYGBgCgoKYGBge3J9CiMjIEZpZy4gMy4zLgpwbG90KHkgfiB4LCBkYXRhID0gQ2h3aXJ1dDIsIHhsYWIgPSAiTWV0YWwgZGlzdGFuY2UiLAogICAgIHlsYWIgPSAiVWx0cmFzb25pYyByZXNwb25zZSIsIHlsaW0gPSBjKDAsIDEwMCkpCmN1cnZlKGV4cEZjdCh4LCBiZXRhMSA9IDAuMSwgYmV0YTIgPSAwLjAxLCBiZXRhMyA9IDEpLCBhZGQgPSBUUlVFLCBsdHkgPSAyLAogICAgICBjb2wgPSAyKQpjdXJ2ZShleHBGY3QoeCwgYmV0YTEgPSAwLjEsIGJldGEyID0gMC4wMSwgYmV0YTMgPSAwLjEpLCBhZGQgPSBUUlVFLCBsdHkgPSAzLAogICAgICBjb2wgPSAzKQpjdXJ2ZShleHBGY3QoeCwgYmV0YTEgPSAwLjEsIGJldGEyID0gMC4wMSwgYmV0YTMgPSAwLjAxKSwgYWRkID0gVFJVRSwgbHR5ID0gNCwKICAgICAgY29sID0gNCkKY3VydmUoZXhwRmN0KHgsIGJldGExID0gMC4yLCBiZXRhMiA9IDAuMDEsIGJldGEzID0gMC4wMSksIGFkZCA9IFRSVUUsIGx0eSA9IDEsCiAgICAgIGNvbCA9IDYpCmBgYAoKCmBgYHtyfQpDaHdpcnV0Mi5tMSA8LSBubHMoCiAgeSB+IGV4cEZjdCh4LCBiZXRhMSwgYmV0YTIsIGJldGEzKSwKICBkYXRhID0gQ2h3aXJ1dDIsCiAgc3RhcnQgPSBsaXN0KGJldGExID0gMC4yLCBiZXRhMiA9IDAuMDEsIGJldGEzID0gMC4wMSkKKQpzdW1tYXJ5KENod2lydXQyLm0xKQpgYGAKCgojIyMgMy4xLjIgU2VhcmNoaW5nIGEgZ3JpZAoKYGBge3J9CiMjIGJhZCBzdGFydGluZyBleGFtcGxlCiMgQ2h3aXJ1dDIubTAgPC0gbmxzKAojICAgeSB+IGV4cEZjdCh4LCBiZXRhMSwgYmV0YTIsIGJldGEzKSwKIyAgIGRhdGEgPSBDaHdpcnV0MiwKIyAgIHN0YXJ0ID0gbGlzdChiZXRhMSA9IDEsIGJldGEyID0gMC4wMSwgYmV0YTMgPSAxKSkKIyBzdW1tYXJ5KENod2lydXQyLm0wKQoKIyMgbmxzTE0gbWF5IHNvbHZlIHRoZSBiYWQgc3RhcnRpbmcKQ2h3aXJ1dDIubTAgPC0gbWlucGFjay5sbTo6bmxzTE0oCiAgeSB+IGV4cEZjdCh4LCBiZXRhMSwgYmV0YTIsIGJldGEzKSwKICBkYXRhID0gQ2h3aXJ1dDIsCiAgc3RhcnQgPSBsaXN0KGJldGExID0gMSwgYmV0YTIgPSAwLjAxLCBiZXRhMyA9IDEpKQpzdW1tYXJ5KENod2lydXQyLm0wKQpgYGAKCgpgYGB7cn0KZ3JpZC5DaHdpcnV0MiA8LSBleHBhbmQuZ3JpZChiZXRhMSA9IHNlcSgwLjEsIDEsIGJ5ID0gMC4xKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBiZXRhMiA9ICAwLjAxLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIGJldGEzID0gc2VxKDAuMSwgMSwgYnkgPSAwLjEpKQpncmlkLkNod2lydXQyCmBgYAoKCmBgYHtyfQpDaHdpcnV0Mi5tMmEgPC0gbmxzMigKICB5IH4gZXhwRmN0KHgsIGJldGExLCBiZXRhMiwgYmV0YTMpLAogIGRhdGEgPSBDaHdpcnV0MiwKICBzdGFydCA9IGdyaWQuQ2h3aXJ1dDIsCiAgYWxnb3JpdGhtID0gImJydXRlLWZvcmNlIgopCkNod2lydXQyLm0yYQpgYGAKCgpgYGB7cn0KbmxzKAogIHkgfiBleHBGY3QoeCwgYmV0YTEsIGJldGEyLCBiZXRhMyksCiAgZGF0YSA9IENod2lydXQyLAogIHN0YXJ0ID0gYXMubGlzdChjb2VmKENod2lydXQyLm0yYSkpKSAgIyBzdGFydGluZyB2YWx1ZSBmcm9tIENod2lydXQyLm0yYQpgYGAKCgojIyAzLjIgVXNpbmcgc2VsZi1zdGFydGVyIGZ1bmN0aW9ucwoKIyMjIDMuMi4xIEJ1aWx0LWluIHNlbGYtc3RhcnRlciBmdW5jdGlvbnMgZm9yIG5scygpCgpgYGB7cn0KTC5taW5vci5tMiA8LSBubHMocmF0ZSB+IFNTbWljbWVuKGNvbmMsIFZtLCBLKSwgZGF0YSA9IEwubWlub3IpCnN1bW1hcnkoTC5taW5vci5tMikKYGBgCgoKIyMjIDMuMi4yIERlZmluaW5nIGEgc2VsZi1zdGFydGVyIGZ1bmN0aW9uIGZvciBubHMoKQoKYGBge3J9CmV4cE1vZGVsIDwtIGZ1bmN0aW9uKHByZWRpY3RvciwgYiwgeTApIHsKICB5MCAqIGV4cChwcmVkaWN0b3IvYikKfQoKZXhwTW9kZWxJbml0IDwtIGZ1bmN0aW9uKG1DYWxsLCBMSFMsIGRhdGEpIHsKICB4eSA8LSBzb3J0ZWRYeURhdGEobUNhbGxbWyJwcmVkaWN0b3IiXV0sIExIUywgZGF0YSkKICBsbUZpdCA8LSBsbShsb2coeHlbLCAieSJdKSB+IHh5WywgIngiXSkKICBjb2VmcyA8LSBjb2VmKGxtRml0KQogIHkwIDwtIGV4cChjb2Vmc1sxXSkKICBiIDwtIDEvY29lZnNbMl0KICB2YWx1ZSA8LSBjKGIsIHkwKQogIG5hbWVzKHZhbHVlKSA8LSBtQ2FsbFtjKCJiIiwgInkwIildCiAgdmFsdWUKfQoKU1NleHAyIDwtIHNlbGZTdGFydChleHBNb2RlbCwgZXhwTW9kZWxJbml0LCBjKCJiIiwgInkwIikpCmBgYAoKCmBgYHtyfQp3aXRoKFJHUmN1cnZlLCBTU2V4cDIoRGF5LCA0LCAwLjIpKQpgYGAKCioqV2FybmluZzoqKiBzZWxmU3RhcnQgY2F1c2VkIGFuIEVSUk9SCgpgYGB7cn0KIyMgRVJST1IKIyBnZXRJbml0aWFsKFJHUiB+IFNTZXhwMihEYXksIGIsIHkwKSwgZGF0YSA9IFJHUmN1cnZlKQpgYGAKCgpgYGB7cn0KI1JHUmN1cnZlLm0xIDwtIG5scyhSR1IgfiBTU2V4cDIoRGF5LCBiLCB5MCksIGRhdGEgPSBSR1JjdXJ2ZSkKI2NvZWYoUkdSY3VydmUubTEpCmBgYAoKTW9yZSBleGFtcGxlcyBvbiBkZWZpbmluZyBzZWxmLXN0YXJ0ZXIgZnVuY3Rpb25zIGZvciBubHMoKSBhcmUgZm91bmQgaW4gV2F0a2lucyBhbmQgVmVuYWJsZXMgKDIwMDYpIG9yIFZlbmFibGVzIGFuZCBSaXBsZXkgKDIwMDJhLCBwcC4gMjE24oCTMjE3KS4KCiMgNCBNb3JlIG9uIG5scygpCgojIyA0LjIgU3VwcGx5aW5nIGdyYWRpZW50IGluZm9ybWF0aW9uCgojIyMgNC4yLjEgTWFudWFsIHN1cHBseQoKYGBge3J9Ck1NZmN0MSA8LSBmdW5jdGlvbihjb25jLCBLLCBWbSkgewogIG51bWVyIDwtIFZtICogY29uYwogIGRlbm9tIDwtIEsgKyBjb25jCiAgbWVhbiA8LSBudW1lci9kZW5vbQogIHBhcnRpYWxLIDwtIC1udW1lci8oZGVub21eMikKICBwYXJ0aWFsVm0gPC0gbWVhbi9WbQogIGF0dHIobWVhbiwgImdyYWRpZW50IikgPC0gY2JpbmQocGFydGlhbEssIHBhcnRpYWxWbSkKICByZXR1cm4obWVhbikKfQpNTWZjdDEoMSwgMiwgMykKYGBgCgoKYGBge3J9CkwubWlub3IubWdyMSA8LSBubHMocmF0ZSB+IE1NZmN0MShjb25jLCBLLCBWbSksIAogICAgICAgICAgICAgICAgICAgIGRhdGEgPSBMLm1pbm9yLCAKICAgICAgICAgICAgICAgICAgICBzdGFydCA9IGxpc3QoSyA9IDEsIFZtID0gMSksCiAgICAgICAgICAgICAgICAgICAgdHJhY2UgPSBUUlVFKQpzdW1tYXJ5KEwubWlub3IubWdyMSkKYGBgCgoKYGBge3J9CiMjIHdpdGhvdXQgZGVyaXZhdGl2ZQpubHMocmF0ZSB+IFZtICogY29uYy8oSyArIGNvbmMpLCAKICAgIGRhdGEgPSBMLm1pbm9yLCBzdGFydCA9IGxpc3QoSyA9IDEsIFZtID0gMSksIAogICAgdHJhY2UgPSBUUlVFKQpgYGAKCgojIyMgNC4yLjIgQXV0b21hdGljIHN1cHBseQoKYGBge3J9Ck1NZmN0MiA8LSBkZXJpdih+IFZtICogY29uYyAvIChLICsgY29uYyksIAogICAgICAgICAgICAgICAgYygiSyIsICJWbSIpLCAKICAgICAgICAgICAgICAgIGZ1bmN0aW9uKGNvbmMsIEssIFZtKSB7fSkKTU1mY3QyCk1NZmN0MigxLCAyLCAzKQpgYGAKCgpgYGB7cn0KTC5taW5vci5tZ3IyIDwtIG5scygKICByYXRlIH4gTU1mY3QyKGNvbmMsICsgSywgVm0pLCAKICBkYXRhID0gTC5taW5vciwgCiAgc3RhcnQgPSBsaXN0KEsgPSAyMCwgVm0gPSAxMjApKQpzdW1tYXJ5KEwubWlub3IubWdyMikKYGBgCgpSZWY6IFZlbmFibGVz7JmAIFJpcGxleSgyMDAyYSwgcHAuIDIxNOKAkzIxNikKCiMjIDQuMyBDb25kaXRpb25hbGx5IGxpbmVhciBwYXJhbWV0ZXJzCgojIyMgNC4zLjEgbmxzKCkgdXNpbmcgdGhlICJwbGluZWFyIiBhbGdvcml0aG0KCmBgYHtyfQpMLm1pbm9yLm0zIDwtIG5scyhyYXRlIH4gY29uYy8oSyArIGNvbmMpLCAKICAgICAgICAgICAgICAgICAgZGF0YSA9IEwubWlub3IsIAogICAgICAgICAgICAgICAgICBhbGdvcml0aG0gPSAicGxpbmVhciIsCiAgICAgICAgICAgICAgICAgIHN0YXJ0ID0gbGlzdChLID0gMjApKQpzdW1tYXJ5KEwubWlub3IubTMpCmBgYAoKCiMjIyA0LjMuMiBBIHBlZGVzdHJpYW4gYXBwcm9hY2gKCmBgYHtyfQojIyBGaWcuIDQuMS4KZGF0YShzZWdyZWcpCnBsb3QoQyB+IFRlbXAsIGRhdGEgPSBzZWdyZWcsIAogICAgIHhsYWIgPSBNZWFuIH4gdGVtcGVyYXR1cmUgfiAoZGVncmVlIH4gRiksCiAgICAgeWxhYiA9ICJFbmVyZ3kgY29uc3VtcHRpb24gKGtXaCkiKQpgYGAKCgpgYGB7cn0KcHJvZlJTUzEgPC0gZnVuY3Rpb24oZ2FtbWEpIHsKICBkZXZpYW5jZShsbShDIH4gcG1heCgwLCBUZW1wIC0gZ2FtbWEpLCBkYXRhID0gc2VncmVnKSkKfQpwcm9mUlNTMiA8LSBWZWN0b3JpemUocHJvZlJTUzEsICJnYW1tYSIpCmBgYAoKCmBgYHtyfQpsIDwtIGxtKCBDIH4gcG1heCgwLCBUZW1wIC0gNDApLCBkYXRhID0gc2VncmVnKQpzdW1tYXJ5KGwpCmBgYAoKCmBgYHtyfQpwbG90KHNlZ3JlZyRUZW1wLCBzZWdyZWckQykKbGluZXMoc2VncmVnJFRlbXAsIHByZWRpY3QobCksIGNvbCA9IDIpCmBgYAoKCmBgYHtyfQojIyBGaWcuIDQuMi4gCnBsb3QocHJvZlJTUzIoVGVtcCkgfiBUZW1wLCAKICAgICBkYXRhID0gc2VncmVnLAogICAgIHR5cGUgPSAibCIsIAogICAgIHhsYWIgPSBleHByZXNzaW9uKGdhbW1hKSwKICAgICB5bGFiID0gIlByb2ZpbGUgUlNTIikKYGBgCgoKYGBge3J9CiMjIOyigOuNlCDrtoDrk5zrn73qsowg66eM65OsCnBsb3QocHJvZlJTUzIoVGVtcCkgfiBUZW1wLCAKICAgICBkYXRhID0gZGF0YS5mcmFtZShUZW1wID0gc2VxKG1pbihzZWdyZWckVGVtcCksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWF4KHNlZ3JlZyRUZW1wKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBsZW5ndGgub3V0ID0gMTAwMSkpLAogICAgIHR5cGUgPSAibCIsIAogICAgIHhsYWIgPSBleHByZXNzaW9uKGdhbW1hKSwKICAgICB5bGFiID0gIlByb2ZpbGUgUlNTIikKYGBgCgoKIyMgNC40IEZpdHRpbmcgbW9kZWxzIHdpdGggc2V2ZXJhbCBwcmVkaWN0b3IgdmFyaWFibGVzCgojIyMgNC40LjEgVHdvLWRpbWVuc2lvbmFsIHByZWRpY3RvcgoKYGBge3J9ClJTY29tcGV0aXRpb24ubTEgPC0gbmxzKAogIGJpb21hc3MgfiBhLygxICsgYiAqICh4ICsgYyAqIHopKSwKICBkYXRhID0gUlNjb21wZXRpdGlvbiwKICBzdGFydCA9IGxpc3QoYSA9IDIwLCBiID0gMSwgYyA9IDEpCikKc3VtbWFyeShSU2NvbXBldGl0aW9uLm0xKQpgYGAKCgpgYGB7cn0KdmlyRGVuc2l0eSA8LSB3aXRoKFJTY29tcGV0aXRpb24sIHggKyBjb2VmKFJTY29tcGV0aXRpb24ubTEpWzNdICogeikKdmlyRGVuVmFsIDwtIHNlcSgwLCBtYXgodmlyRGVuc2l0eSksIGxlbmd0aC5vdXQgPSAxMDApCmJpb21hc3NWYWwgPC0gcHJlZGljdChSU2NvbXBldGl0aW9uLm0xLCBkYXRhLmZyYW1lKHggPSB2aXJEZW5WYWwsIHogPSAwKSkKYGBgCgoKYGBge3J9CiMjIEZpZy4gNC4zLiAKcGxvdChiaW9tYXNzVmFsIH4gdmlyRGVuVmFsLCB0eXBlID0gImwiLCAKICAgICB5bGFiID0gIkJpb21hc3Mgb2Ygc2Vuc2l0aXZlIGJpb3R5cGUgKGcvcGxhbnQpIiwKICAgICB4bGFiID0gVmlydHVhbCB+IGRlbnNpdHkgfiAocGxhbnRzL21eMikpCndpdGgoUlNjb21wZXRpdGlvbiwgcG9pbnRzKGJpb21hc3MgfiB2aXJEZW5zaXR5KSkKYGBgCgoKIyMjIDQuNC4yIEdlbmVyYWwgbGVhc3Qtc3F1YXJlcyBtaW5pbWlzYXRpb24KCmBgYHtyfQojIyBGaWcuIDQuNC4gCnBhcihwdHkgPSAicyIpCnBsb3QoUSB+IEksIGRhdGEgPSBJUXNpZykKdGhldGEgPC0gMDozNjAgKiAocGkvMTgwKQpsaW5lcyhjb3ModGhldGEpLCBzaW4odGhldGEpKQpgYGAKCgpgYGB7cn0KSVFzaWcubTEgPC0gbmxzKAogIH4gKChJIC0gSTApXjIgLSAyICogZ2FtbWEgKiBzaW4ocGhpKSAqIChJIC0gSTApICogKFEgLSBRMCkgKyAKICAgICAgIGdhbW1hICogZ2FtbWEgKiAoUSAtIFEwKV4yKSAtIChyaG8gKiBnYW1tYSAqIGNvcyhwaGkpKV4yLCAKICBkYXRhID0gSVFzaWcsIAogIHN0YXJ0ID0gbGlzdChJMCA9IC0wLjAwNSwgZ2FtbWEgPSAxLCBwaGkgPSAtMC4wMDUsIFEwID0gLTAuMDA1LCByaG8gPSAxKSkKc3VtbWFyeShJUXNpZy5tMSkKYGBgCgoKIyA1IE1vZGVsIERpYWdub3N0aWNzCgojIyA1LjIgQ2hlY2tpbmcgdGhlIG1lYW4gc3RydWN0dXJlCgojIyMgNS4yLjEgUGxvdCBvZiB0aGUgZml0dGVkIHJlZ3Jlc3Npb24gY3VydmUKCmBgYHtyfQojIyBGaWcuIDUuMS4KcGxvdChwIH4gVCwgZGF0YSA9IHZhcENPLCBsb2cgPSAieSIsCiAgICAgeGxhYiA9ICJUZW1wZXJhdHVyZSAoSykiLAogICAgIHlsYWIgPSAiUHJlc3N1cmUgKFBhKSIpCmBgYAoKCmBgYHtyfQp2YXBDTy5tMSA8LSBubHMoCiAgbG9nKHApIH4gQSAtIEIvKEMgKyBUKSwgCiAgZGF0YSA9IHZhcENPLAogIHN0YXJ0ID0gbGlzdChBID0gMTAsIEIgPSAxMDAsIEMgPSAtMTApKQpzdW1tYXJ5KHZhcENPLm0xKQpgYGAKCgpgYGB7cn0KIyMgRmlnLiA1LjIuCnBsb3QocCB+IFQsIGRhdGEgPSB2YXBDTywgbG9nID0gInkiLAogICAgIHhsYWIgPSAiVGVtcGVyYXR1cmUgKEspIiwKICAgICB5bGFiID0gIlByZXNzdXJlIChQYSkiKQpsaW5lcyh2YXBDTyRULCBleHAoZml0dGVkKHZhcENPLm0xKSksIGNvbCA9IDIpCmBgYAoKCmBgYHtyfQojIyBGaWcuIDUuMy4gCnBsb3Qod2VpZ2h0IH4gY29uYywgZGF0YSA9IGxldHR1Y2UsCiAgICAgeGxhYiA9ICJDb25jZW50cmF0aW9uIChtZy9sKSIsCiAgICAgeWxhYiA9ICJCaW9tYXNzIChnKSIsIGxvZyA9ICJ4IikKYGBgCgoKIyMjIDUuMi4yIFJlc2lkdWFsIHBsb3RzCgpgYGB7cn0KIyMgRmlnLiA1LjQuCnBsb3QoZml0dGVkKHZhcENPLm0xKSwgcmVzaWR1YWxzKHZhcENPLm0xKSwKICAgICB4bGFiID0gIkZpdHRlZCBWYWx1ZXMiLCB5bGFiID0gIlJlc2lkdWFscyIpCmFibGluZShhID0gMCwgYiA9IDApCmBgYAoKCmBgYHtyfQojIyBGaWcuIDUuNS4gCnBsb3Qocm9vdGwgfiBjb25jLCBkYXRhID0gcnllZ3Jhc3MsCiAgICAgeGxhYiA9ICJDb25jZW50cmF0aW9uIChtTSkiLAogICAgIHlsYWIgPSAiUm9vdCBsZW5ndGggKGNtKSIpCmBgYAoKCmBgYHtyfQpyeWVncmFzcy5tMSA8LSBsbShyb290bCB+IGFzLmZhY3Rvcihjb25jKSwgZGF0YSA9IHJ5ZWdyYXNzKQpzdW1tYXJ5KHJ5ZWdyYXNzLm0xKQpgYGAKCgpgYGB7cn0KcnllZ3Jhc3MubTIgPC0gbmxzKAogIHJvb3RsIH4gYyArIChkIC0gYykvKDEgKyBleHAoYiAqICsobG9nKGNvbmMpIC0gbG9nKGUpKSkpLAogIGRhdGEgPSByeWVncmFzcywKICBzdGFydCA9IGxpc3QoYiA9IDEsIGMgPSAwLjYsIGQgPSA4LCBlID0gMykKKQpzdW1tYXJ5KHJ5ZWdyYXNzLm0yKQpgYGAKCgpgYGB7cn0KYW5vdmEocnllZ3Jhc3MubTIsIHJ5ZWdyYXNzLm0xKQpgYGAKCgpgYGB7cn0KUSA8LSAtMiAqIChsb2dMaWsocnllZ3Jhc3MubTIpIC0gbG9nTGlrKHJ5ZWdyYXNzLm0xKSkKZGYuUSA8LSBkZi5yZXNpZHVhbChyeWVncmFzcy5tMikgLSBkZi5yZXNpZHVhbChyeWVncmFzcy5tMSkKMSAtIHBjaGlzcShRLCBkZi5RKQpgYGAKCgojIyA1LjMgVmFyaWFuY2UgaG9tb2dlbmVpdHkKCiMjIyA1LjMuMiBMZXZlbmXigJlzIHRlc3QKCmBgYHtyfQp3aXRoKHJ5ZWdyYXNzLCBsZXZlbmVUZXN0KHJvb3RsLCBhcy5mYWN0b3IoY29uYykpKQpgYGAKCgpgYGB7cn0KIyMgRmlnLiA1LjYuIApwbG90KGZpdHRlZCh2YXBDTy5tMSksIGFicyhyZXNpZHVhbHModmFwQ08ubTEpKSwKICAgICB4bGFiID0gIkZpdHRlZCB2YWx1ZXMiLCB5bGFiID0gIkFic29sdXRlIHJlc2lkdWFscyIpCmBgYAoKCiMjIDUuNCBOb3JtYWwgZGlzdHJpYnV0aW9uCgpgYGB7cn0KIyMgRmlnLiA1LjcuIApwbG90KHZhcENPLm0xKQpgYGAKCgojIyMgNS40LjEgUVEgcGxvdAoKYGBge3J9CiMjIEZpZy4gNS44LiAKc3RhbmRhcmRSZXMgPC0gcmVzaWR1YWxzKHJ5ZWdyYXNzLm0yKS9zdW1tYXJ5KHJ5ZWdyYXNzLm0yKSRzaWdtYQoKcGFyKHB0eSA9ICJzIikKcXFub3JtKHN0YW5kYXJkUmVzLCBtYWluID0gIiIpCmFibGluZShhID0gMCwgYiA9IDEpCmBgYAoKCiMjIyA1LjQuMiBTaGFwaXJvLVdpbGsgdGVzdAoKYGBge3J9CnNoYXBpcm8udGVzdChzdGFuZGFyZFJlcykKYGBgCgoKIyMgNS41IEluZGVwZW5kZW5jZQoKYGBge3J9CiMjIEZpZy4gNS45LiAKcGxvdChyZXNpZHVhbHModmFwQ08ubTEpLCBjKHJlc2lkdWFscyh2YXBDTy5tMSlbLTFdLCBOQSksCiAgICAgeGxhYiA9ICJSZXNpZHVhbHMiLCB5bGFiID0gIkxhZ2dlZCByZXNpZHVhbHMiKQpgYGAKCgojIDYgUmVtZWRpZXMgZm9yIE1vZGVsIFZpb2xhdGlvbnMKCiMjIDYuMSBWYXJpYW5jZSBtb2RlbGxpbmcKCiMjIyA2LjEuMSBQb3dlci1vZi10aGUtbWVhbiB2YXJpYW5jZSBtb2RlbAoKYGBge3J9CiMjIEZpZy4gNi4xLiAKcGxvdChSR1IgfiBEYXksIGRhdGEgPSBSR1JjdXJ2ZSwKICAgICB4bGFiID0gIlRpbWUgKGRheXMpIiwgeWxhYiA9ICJSZWxhdGl2ZSBncm93dGggcmF0ZSAoJSkiKQpgYGAKCgpgYGB7cn0KeTBfc3RhcnQgPSBtZWFuKFJHUmN1cnZlW1JHUmN1cnZlJERheSA9PSAwLCAiUkdSIl0pCnkxID0gbWVhbihSR1JjdXJ2ZVtSR1JjdXJ2ZSREYXkgPT0gMSwgIlJHUiJdKQpiX3N0YXJ0ID0gMSAvIGxvZyh5MS95MF9zdGFydCkKClJHUmN1cnZlLm0yIDwtIGdubHMoCiAgUkdSIH4gIHkwICogZXhwKERheS9iKSwgCiAgZGF0YSA9IFJHUmN1cnZlLAogIHN0YXJ0ID0gbGlzdCh5MCA9IHkwX3N0YXJ0LCBiID0gYl9zdGFydCksCiAgd2VpZ2h0cyA9IHZhclBvd2VyKCkKKQpzdW1tYXJ5KFJHUmN1cnZlLm0yKQpgYGAKCgpgYGB7cn0KUkdSY3VydmUubTJnIDwtIGdubHMoCiAgUkdSIH4gIHkwICogZXhwKERheS9iKSwgCiAgZGF0YSA9IFJHUmN1cnZlLAogIHN0YXJ0ID0gbGlzdCh5MCA9IHkwX3N0YXJ0LCBiID0gYl9zdGFydCkKICAjIHdlaWdodHMgPSB2YXJQb3dlcigpICMgd2l0aG91dCB3ZWlnaHRzCikKc3VtbWFyeShSR1JjdXJ2ZS5tMmcgKQpgYGAKCgpgYGB7cn0KQUlDKFJHUmN1cnZlLm0yLCBSR1JjdXJ2ZS5tMmcgKQpgYGAKCgojIyA2LjIgVHJhbnNmb3JtYXRpb25zCgojIyMgNi4yLjEgVHJhbnNmb3JtLWJvdGgtc2lkZXMgYXBwcm9hY2gKCmBgYHtyfQojIyBGaWcuIDYuMi4gCnBsb3QocmVjcnVpdHMgfiBzcGF3bmVycywgCiAgICAgZGF0YSA9IHNvY2tleWVbLTEyLF0sCiAgICAgeGxhYiA9ICJOdW1iZXIgb2Ygc3Bhd25lcnMgKHRob3VzYW5kcykiLAogICAgIHlsYWIgPSAiTnVtYmVyIG9mIHJlY3J1aXRzICh0aG91c2FuZHMpIikKYGBgCgoK7J2864uoIG5sc+ulvCDsnbTsmqntlZjsl6wgZml0dGluZ+ydhCDtlZwg7ZuECgpgYGB7cn0Kc29ja2V5ZS5tMSA8LSBubHMoCiAgcmVjcnVpdHMgfiBiZXRhMSAqIHNwYXduZXJzICogZXhwKC1iZXRhMiAqIHNwYXduZXJzKSwKICBkYXRhID0gc29ja2V5ZVstMTIsIF0sCiAgc3RhcnQgPSBsaXN0KGJldGExID0gMiwgYmV0YTIgPSAwLjAwMSkKKQpzdW1tYXJ5KHNvY2tleWUubTEpCmBgYAoKCmBgYHtyfQojIyBGaWcuIDYuMy4gCnBsb3QoZml0dGVkKHNvY2tleWUubTEpLCBhYnMocmVzaWR1YWxzKHNvY2tleWUubTEpKSwKICAgICB4bGFiID0gIkZpdHRlZCB2YWx1ZXMiLCB5bGFiID0gIkFic29sdXRlIHJlc2lkdWFscyIpCmBgYAoKCiMjIyA2LjIuMiBGaW5kaW5nIGFuIGFwcHJvcHJpYXRlIHRyYW5zZm9ybWF0aW9uCgpib3hjb3gubmxz66W8IOydtOyaqe2VmOyXrCDslpHsqr3snYQg67OA7ZmU7Iuc7YKk66mwIO2UvO2MheydhCDtlZzri6QuCgpgYGB7cn0KIyMgRmlnLiA2LjQuCnNvY2tleWUubTIgPC0gYm94Y294Lm5scyhzb2NrZXllLm0xKQpgYGAKCgpgYGB7cn0KYmNTdW1tYXJ5KHNvY2tleWUubTIpCmBgYAoKCmBgYHtyfQpjb2VmKHN1bW1hcnkoc29ja2V5ZS5tMSkpCmNvZWYoc3VtbWFyeShzb2NrZXllLm0yKSkKYGBgCgoKIyMgNi4zIFNhbmR3aWNoIGVzdGltYXRvcnMKCuydtCDrgrTsmqnsnYAg7KCV7ZmV7Z6IIOydtO2VtO2VmOyngOuKlCDrqrvtlojsnYwuCgpgYGB7cn0KdmNvdihzb2NrZXllLm0xKQpgYGAKCgpgYGB7cn0Kc2FuZHdpY2goc29ja2V5ZS5tMSkKYGBgCgoKYGBge3J9CmNvZWZ0ZXN0KHNvY2tleWUubTEpCmBgYAoKCmBgYHtyfQpjb2VmdGVzdChzb2NrZXllLm0xLCB2Y292ID0gc2FuZHdpY2gpCmBgYAoKCiMjIDYuNCBXZWlnaHRpbmcKCiMjIyA2LjQuMSBEZWNsaW5lIGluIG5pdHJvZ2VuIGNvbnRlbnQgaW4gc29pbAoKYGBge3J9CmV4cDEKYGBgCgoKYGBge3J9CiMjIEZpZy4gNi41LiAKcGxvdChOcmVtYWluaW5nIH4gdGltZSwgZGF0YSA9IGV4cDEsCiAgICAgeGxhYiA9ICJUaW1lICh5ZWFycykiLCB5bGFiID0gIk5pdHJvZ2UgY29udGVudCAoJSkiKQpgYGAKCgpgYGB7cn0KZXhwMS5tMSA8LSBubHMoTnJlbWFpbmluZyB+IFNTYmlleHAodGltZSwgYTEsIGEyLCBiMSwgYjIpLCBkYXRhID0gZXhwMSkKZXhwMS5tMiA8LSBubHMoTnJlbWFpbmluZyB+IFNTYmlleHAodGltZSwgYTEsIGEyLCBiMSwgYjIpLCBkYXRhID0gZXhwMSwKICAgICAgICAgICAgICAgd2VpZ2h0cyA9IG5vcmVwLyhzdGRldiAqIHN0ZGV2KSkKYGBgCgoKYGBge3J9CndlaWdodHMoZXhwMS5tMikKYGBgCgoKYGBge3J9CmNvZWYoc3VtbWFyeShleHAxLm0xKSkKY29lZihzdW1tYXJ5KGV4cDEubTIpKQpgYGAKCgpgYGB7cn0KZXhwMgpgYGAKCgpgYGB7cn0KZXhwMi5tMSA8LSBubHMoTnJlbWFpbmluZyB+IFNTYmlleHAodGltZSwgYTEsIGEyLCBiMSwgYjIpLCBkYXRhID0gZXhwMikKZXhwMi5tMiA8LSBubHMoTnJlbWFpbmluZyB+IFNTYmlleHAodGltZSwgYTEsIGEyLCBiMSwgYjIpLCBkYXRhID0gZXhwMiwKICAgICAgICAgICAgICAgd2VpZ2h0cyA9IG5vcmVwLyhzdGRldl4yKSkKYGBgCgoKYGBge3J9CiMjIEZpZy4gNi42LiAKcGxvdChOcmVtYWluaW5nIH4gdGltZSwgZGF0YSA9IGV4cDIsCiAgICAgeGxhYiA9ICJUaW1lICh5ZWFycykiLCB5bGFiID0gIk5pdHJvZ2VuIGNvbnRlbnQgKCUpIikKdGltZVZhbCA8LSB3aXRoKGV4cDIsIHNlcShtaW4odGltZSksIG1heCh0aW1lKSwgbGVuZ3RoLm91dCA9IDEwMSkpCmxpbmVzKHRpbWVWYWwsIHByZWRpY3QoZXhwMi5tMSwgbmV3ZGF0YSA9IGRhdGEuZnJhbWUodGltZSA9IHRpbWVWYWwpKSwgCiAgICAgIGx0eSA9IDIsIGNvbCA9IDIpCmxpbmVzKHRpbWVWYWwsIHByZWRpY3QoZXhwMi5tMiwgbmV3ZGF0YSA9IGRhdGEuZnJhbWUodGltZSA9IHRpbWVWYWwpKSwgCiAgICAgIGx0eSA9IDMsIGNvbCA9IDMpCmBgYAoKCmBgYHtyfQpjb2VmKHN1bW1hcnkoZXhwMi5tMSkpCmNvZWYoc3VtbWFyeShleHAyLm0yKSkKYGBgCgoKYGBge3J9CnN1bW1hcnkoZXhwMS5tMSkKc3VtbWFyeShleHAxLm0yKQpgYGAKCgojIDcgVW5jZXJ0YWludHksIEh5cG90aGVzaXMgVGVzdGluZywgYW5kIE1vZGVsIFNlbGVjdGlvbgoKIyMgNy4xIFByb2ZpbGUgbGlrZWxpaG9vZAoKYGBge3J9CkwubWlub3IubTEgPC0gdXBkYXRlKEwubWlub3IubTEsIHRyYWNlID0gRkFMU0UpCkwubWlub3IubTFwcm8gPC0gcHJvZmlsZShMLm1pbm9yLm0xKQpgYGAKCgpgYGB7cn0KIyMgRmlnLiA3LjEuIApwbG90KEwubWlub3IubTFwcm8pCmBgYAoKCmBgYHtyfQojIyBGaWcuIDcuMi4gCnBsb3QoTC5taW5vci5tMXBybywgYWJzVmFsID0gRkFMU0UpCmBgYAoKCmBgYHtyfQpjb25maW50KEwubWlub3IubTEpCmBgYAoKCiMjIDcuMiBCb290c3RyYXAKCmBgYHtyfQpzZXQuc2VlZCgyMjIpCkwubWlub3IubTFib290IDwtIG5sc0Jvb3QoTC5taW5vci5tMSkKc3VtbWFyeShMLm1pbm9yLm0xYm9vdCkKYGBgCgoKYGBge3J9CiMjIEZpZy4gNy4zLiAKcGFyKG1mcm93ID0gYygxLCAyKSkKcXFub3JtKEwubWlub3IubTFib290JGNvZWZib290WywgMV0sIG1haW4gPSAiSyIpCnFxbm9ybShMLm1pbm9yLm0xYm9vdCRjb2VmYm9vdFssIDJdLCBtYWluID0gIlZtIikKYGBgCgoKIyMgNy4zIFdhbGQgY29uZmlkZW5jZSBpbnRlcnZhbHMKCmBgYHtyfQpjb25maW50MihMLm1pbm9yLm0xKQpgYGAKCgpgYGB7cn0KY29uZmludDIoTC5taW5vci5tMSwgbGV2ZWwgPSAwLjk5KQpgYGAKCgojIyA3LjQgRXN0aW1hdGluZyBkZXJpdmVkIHBhcmFtZXRlcnMKCmBgYHtyfQpkZWx0YU1ldGhvZChMLm1pbm9yLm0xLCAiVm0vKDQqSykiKQpgYGAKCgojIyA3LjUgTmVzdGVkIG1vZGVscwoKYGBge3J9CnNlY2Fsb25pYwpgYGAKCgpgYGB7cn0KIyMgRmlnLiA3LjQuCnBsb3Qocm9vdGwgfiBkb3NlLCBkYXRhID0gc2VjYWxvbmljLAogICAgIHhsYWIgPSAiRG9zZSAobU0pIiwgeWxhYiA9ICJSb290IGxlbmd0aCAoY20pIikKYGBgCgoKIyMjIDcuNS4xIFVzaW5nIHQtdGVzdHMKCmBgYHtyfQpzZWNhbG9uaWMubTEgPC0gbmxzKHJvb3RsIH4gU1NmcGwoZG9zZSwgYSwgYiwgYywgZCksIGRhdGEgPSBzZWNhbG9uaWMpCnN1bW1hcnkoc2VjYWxvbmljLm0xKQpgYGAKCgojIyMgNy41LjIgVXNpbmcgRi10ZXN0cwoKYGBge3J9CnNlY2Fsb25pYy5tMiA8LSBubHMocm9vdGwgfiBTU2xvZ2lzKGRvc2UsIGEsIGMsIGQpLCBkYXRhID0gc2VjYWxvbmljKQpzdW1tYXJ5KHNlY2Fsb25pYy5tMikKYGBgCgoKYGBge3J9CmFub3ZhKHNlY2Fsb25pYy5tMiwgc2VjYWxvbmljLm0xKQpgYGAKCgojIyA3LjYgTm9uLW5lc3RlZCBtb2RlbHMKCmBgYHtyfQpNLm1lcmx1Y2NpdXMKcGxvdChudW0uZmlzaCB+IHNwYXduLmJpb21hc3MsIGRhdGEgPSBNLm1lcmx1Y2NpdXMpCmBgYAoKNOqwgOyngCDsobDquIjslKkg64uk66W4IOuqqOuNuOydhCDruYTqtZAKCmBgYHtyfQpNLm1lcmx1Y2NpdXMuYmggPC0gbmxzKAogIG51bS5maXNoIH4gc3Bhd24uYmlvbWFzcyAqIGFscGhhLygxICsgc3Bhd24uYmlvbWFzcy9rKSwKICBkYXRhID0gTS5tZXJsdWNjaXVzLCBzdGFydCA9IGxpc3QoYWxwaGEgPSA1LCBrID0gNTApCikKCk0ubWVybHVjY2l1cy5kZSA8LSBubHMoCiAgbnVtLmZpc2ggfiBzcGF3bi5iaW9tYXNzICogYWxwaGEgKiAoMSAtIGMgKiBzcGF3bi5iaW9tYXNzL2spXigxL2MpLAogIGRhdGEgPSBNLm1lcmx1Y2NpdXMsIHN0YXJ0ID0gbGlzdChhbHBoYSA9IDQuNCwgayA9IDEwNiwgYyA9IDAuODYpCikKCk0ubWVybHVjY2l1cy5yaSA8LSBubHMoCiAgbnVtLmZpc2ggfiBzcGF3bi5iaW9tYXNzICogYWxwaGEgKiBleHAoLXNwYXduLmJpb21hc3MvayksCiAgZGF0YSA9IE0ubWVybHVjY2l1cywgc3RhcnQgPSBsaXN0KGFscGhhID0gNSwgayA9IDUwKQopCgpNLm1lcmx1Y2NpdXMuc2ggPC0gbmxzKAogIG51bS5maXNoIH4gc3Bhd24uYmlvbWFzcyAqIGFscGhhLygxICsgKHNwYXduLmJpb21hc3MvayleYyksCiAgZGF0YSA9IE0ubWVybHVjY2l1cywgc3RhcnQgPSBsaXN0KGFscGhhID0gMy44NywgayA9IDYxLjcyLCBjID0gMi4yNSksCiAgY29udHJvbCA9IG5scy5jb250cm9sKG1heGl0ZXIgPSAxMDApCikKYGBgCgpSZXNpZHVhbCBzdGFuZGFyZCBlcnJvcuulvCDqs4TsgrDtlZjrqbQsCgpgYGB7cn0Kc3VtbWFyeShNLm1lcmx1Y2NpdXMuYmgpJHNpZ21hCnN1bW1hcnkoTS5tZXJsdWNjaXVzLmRlKSRzaWdtYQpzdW1tYXJ5KE0ubWVybHVjY2l1cy5yaSkkc2lnbWEKc3VtbWFyeShNLm1lcmx1Y2NpdXMuc2gpJHNpZ21hCmBgYAoKQUlD66W8IOu5hOq1kO2VtOuztOuptCwKCmBgYHtyfQpBSUMoTS5tZXJsdWNjaXVzLmJoKQpBSUMoTS5tZXJsdWNjaXVzLmRlKQpBSUMoTS5tZXJsdWNjaXVzLnJpKQpBSUMoTS5tZXJsdWNjaXVzLnNoKQpgYGAKCkJ1cm5oYW0gYW5kIEFuZGVyc29uICgyMDAyLCBwcC4gNzDigJM3MinsnYAg7ZWcIOuqqOuNuOydtCDri6Trpbgg66qo64247JeQIOu5hO2VmOyXrCDsoIjrjIDsoIHsnLzroZwg7Jqw7JyE7JeQIOyeiOq4sCDsnITtlbTshJzripQgQUlD7J2YIOywqOydtOqwgCAxMCDsnbTsg4HsnZgg7LCo7J206rCAIO2VhOyalO2VmOuLpOuKlCDqsr3tl5jsoIEg67KV7LmZ7J2EIOygnOyLnO2VmOyYgOuLpC4g65Sw65287IScIDQg66qo6424IOyCrOydtOyXkOyEnCDtgbAg7LCo7J2064qUIOyXhuuLpC4gCgoKIyA4IEdyb3VwZWQgRGF0YQoKIyMgOC4xIEZpdHRpbmcgZ3JvdXBlZCBkYXRhIG1vZGVscwoKYGBge3J9ClB1cm9teWNpbgpgYGAKCgpgYGB7cn0KIyMgRmlnLiA4LjEuCnh5cGxvdChyYXRlIH4gY29uYyB8IHN0YXRlLCBkYXRhID0gUHVyb215Y2luLAogICAgICAgeGxhYiA9ICJTdWJzdHJhdGUgY29uY2VudHJhdGlvbiAocHBtKSIsCiAgICAgICB5bGFiID0gIlJlYWN0aW9uIHJhdGVzXG4oY291bnRzL21pbi9taW4pIikKYGBgCgoKIyMjIDguMS4xIFVzaW5nIG5scygpCgpgYGB7cn0KUHVyb215Y2luLm0xIDwtIG5scygKICByYXRlIH4gVm1bc3RhdGVdICogY29uYy8oS1tzdGF0ZV0gKyBjb25jKSwKICBkYXRhID0gUHVyb215Y2luLAogIHN0YXJ0ID0gbGlzdChLID0gYygwLjEsIDAuMSksIFZtID0gYygyMDAsIDIwMCkpCikKc3VtbWFyeShQdXJvbXljaW4ubTEpCmBgYAoKCiMjIyA4LjEuMiBVc2luZyBnbmxzKCkKCmBgYHtyfQpQdXJvbXljaW4ubTIgPC0gZ25scygKICByYXRlIH4gVm0gKiBjb25jLyhLICsgY29uYyksCiAgZGF0YSA9IFB1cm9teWNpbiwKICBzdGFydCA9IGxpc3QoVm0gPSBjKDIwMCwgMjAwKSwgSyA9IGMoMC4xLCAwLjEpKSwKICBwYXJhbXMgPSBsaXN0KFZtIH4gc3RhdGUgLSAxLCBLIH4gc3RhdGUgLSAxKQopCnN1bW1hcnkoUHVyb215Y2luLm0yKQpgYGAKCiMjIyA4LjEuMyBVc2luZyBubHNMaXN0KCkKCmBgYHtyfQpQdXJvbXljaW4ubTMgPC0gbmxzTGlzdCgKICByYXRlIH4gU1NtaWNtZW4oY29uYywgVm0sIEspIHwgc3RhdGUsCiAgZGF0YSA9IFB1cm9teWNpbgopCnN1bW1hcnkoUHVyb215Y2luLm0zKQpgYGAKCgpgYGB7cn0KUHVyb215Y2luMiA8LSBncm91cGVkRGF0YShyYXRlIH4gY29uYyB8IHN0YXRlLCBkYXRhID0gUHVyb215Y2luKQpQdXJvbXljaW4ubTQgPC0gbmxzTGlzdCgKICByYXRlIH4gU1NtaWNtZW4oY29uYywgYSwgYiksIGRhdGEgPSBQdXJvbXljaW4yCikKc3VtbWFyeShQdXJvbXljaW4ubTQpCmBgYAoKCiMjIDguMiBNb2RlbCByZWR1Y3Rpb24gYW5kIHBhcmFtZXRlciBtb2RlbHMKCiMjIyA4LjIuMSBDb21wYXJpc29uIG9mIGVudGlyZSBncm91cHMKCmBgYHtyfQpQdXJvbXljaW4ubTUgPC0gbmxzKAogIHJhdGUgfiBWbSAqIGNvbmMvKEsgKyBjb25jKSwgZGF0YSA9IFB1cm9teWNpbiwKICBzdGFydCA9IGxpc3QoSyA9IDAuMSwgVm0gPSAyMDApKQpzdW1tYXJ5KFB1cm9teWNpbi5tNSkKYGBgCgoKYGBge3J9CmFub3ZhKFB1cm9teWNpbi5tNSwgUHVyb215Y2luLm0xKQpgYGAKCgpgYGB7cn0KQUlDKFB1cm9teWNpbi5tNSwgUHVyb215Y2luLm0xKQpgYGAKCiMjIyA4LjIuMiBDb21wYXJpc29uIG9mIHNwZWNpZmljIHBhcmFtZXRlcnMKCmBgYHtyfQpQdXJvbXljaW4ubTYgPC0gbmxzKAogIHJhdGUgfiBWbSAqIGNvbmMvKEtbc3RhdGVdICsgY29uYyksIAogIGRhdGEgPSBQdXJvbXljaW4sIAogIHN0YXJ0ID0gbGlzdChLID0gYygwLjEsIDAuMSksIFZtID0gMjAwKSkKc3VtbWFyeShQdXJvbXljaW4ubTYpCmBgYAoKCmBgYHtyfQphbm92YShQdXJvbXljaW4ubTYsIFB1cm9teWNpbi5tMSkKYGBgCgoKYGBge3J9ClB1cm9teWNpbi5tNyA8LSBubHMoCiAgcmF0ZSB+IFZtW3N0YXRlXSAqIGNvbmMvKEsgKyBjb25jKSwgCiAgZGF0YSA9IFB1cm9teWNpbiwgc3RhcnQgPSBsaXN0KEsgPSAwLjEsIFZtID0gYygyMDAsIDIwMCkpCikKc3VtbWFyeShQdXJvbXljaW4ubTcpCmBgYAoKCmBgYHtyfQphbm92YShQdXJvbXljaW4ubTcsIFB1cm9teWNpbi5tMSkKYGBgCgoKQUlD66GcIOu5hOq1kO2VtCDrs7TslYTrj4QgUHVyb215Y2luLm037J20IOuNlOyasSDsoovri6QuIAoKYGBge3J9CkFJQyhQdXJvbXljaW4ubTEsIFB1cm9teWNpbi5tNSwgUHVyb215Y2luLm02LCBQdXJvbXljaW4ubTcpCmBgYAoKIyMgOC4zIENvbW1vbiBjb250cm9sCgpgYGB7cn0KRy5hcGFyaW5lCmBgYAoKCmBgYHtyfQojIyBGaWcuIDguMi4gCnh5cGxvdChkcnltYXR0ZXIgfiBkb3NlIHwgYXMuZmFjdG9yKHRyZWF0bWVudCksCiAgICAgICBkYXRhID0gRy5hcGFyaW5lLAogICAgICAgeGxhYiA9ICJEb3NlIChnL2hhKSIsIHlsYWIgPSAiRHJ5IG1hdHRlciAobWcvcG90KSIpCmBgYAoK7J2066+4IOuNsOydtO2EsOyXkOyEnCB0cmVhdG1lbnQgMOqwgCDsl4bslrTsoYzsnYwuIOyVhOuniOuPhCDrjbDsnbTthLDrpbwg7JeF642w7J207Yq4IO2VnCDqsoPsnLzroZwg7LaU7Lih65CoLiAKCmBgYHtyfQpHLmFwYXJpbmUkdHJlYXRtZW50MiA8LSBmYWN0b3IoRy5hcGFyaW5lJHRyZWF0bWVudCkKI2xldmVscyhHLmFwYXJpbmUkdHJlYXRtZW50MikgPC0gYygiMCIsICIxIiwgIjIiKQpHLmFwYXJpbmUKYGBgCgoKYGBge3J9CkcuYXBhcmluZS5tMSA8LSBubHMoCiAgZHJ5bWF0dGVyIH4gY1t0cmVhdG1lbnQyXSArIChkIC0gY1t0cmVhdG1lbnQyXSkvKDEgKyBleHAoYlt0cmVhdG1lbnQyXSAqIChsb2coZG9zZSkgLSBsb2coZVt0cmVhdG1lbnQyXSkpKSksCiAgZGF0YSA9IEcuYXBhcmluZSwgCiAgc3RhcnQgPSBsaXN0KGIgPSBjKDIsIDIpLCBjID0gYyg1MDAsIDEwMCksIGQgPSAxMDAwLCBlID0gYyg1MCwgMTAwKSkpCnN1bW1hcnkoRy5hcGFyaW5lLm0xKQpgYGAKCgojIyA4LjQgUHJlZGljdGlvbgoKYGBge3J9CmNvbmNWYWx1ZXMgPC0gd2l0aChQdXJvbXljaW4sIHNlcShtaW4oY29uYyksIG1heChjb25jLCBsZW5ndGgub3V0ID0gMTEpKSkKY29uY1ZhbHVlcwpgYGAKCgpgYGB7cn0Kc3RhdGVWYWwxIDwtIGxldmVscyhQdXJvbXljaW4kc3RhdGUpCnN0YXRlVmFsMQpgYGAKCgpgYGB7cn0KY3NWYWx1ZXMxIDwtIGV4cGFuZC5ncmlkKGNvbmMgPSBjb25jVmFsdWVzLCBzdGF0ZSA9IHN0YXRlVmFsMSkKY3NWYWx1ZXMxCmBgYAoKCmBgYHtyfQpwcmVkaWN0KFB1cm9teWNpbi5tNywgbmV3ZGF0YSA9IGNzVmFsdWVzMSkKYGBgCgoKYGBge3J9CnN0YXRlVmFsMiA8LSBmYWN0b3IoInVudHJlYXRlZCIsIGxldmVscyA9IGMoInRyZWF0ZWQiLCAidW50cmVhdGVkIikpCnN0YXRlVmFsMgpgYGAKCgpgYGB7cn0KY3NWYWx1ZXMyIDwtIGRhdGEuZnJhbWUoY29uYyA9IGNvbmNWYWx1ZXMsIHN0YXRlID0gc3RhdGVWYWwyKQpjc1ZhbHVlczIKYGBgCgoKYGBge3J9CnByZWRpY3QoUHVyb215Y2luLm0xLCBuZXdkYXRhID0gY3NWYWx1ZXMyKQpgYGAKCgojIyA4LjUgTm9ubGluZWFyIG1peGVkIG1vZGVscwoKYGBge3J9CnZpbmNsb3pvbGluCmBgYAoKCmBgYHtyfQojIyBGaWcuIDguMy4KeHlwbG90KGVmZmVjdCB+IGNvbmMgfCBleHBlciwgZGF0YSA9IHZpbmNsb3pvbGluLAogICAgICAgeGxhYiA9IGV4cHJlc3Npb24ocGFzdGUoIkNvbmNlbnRyYXRpb24gKCIsIG11LCAiTSkiKSksCiAgICAgICB5bGFiID0gIkx1bWluZXNjZW5jZSAoTFUiKQpgYGAKCgpgYGB7cn0KTEwzLmZvcm11bGEgPC0gZWZmZWN0IH4gZC8oMSArIGV4cChiICogKGxvZyhjb25jKSAtIGxvZyhlKSkpKQoKdmluY2xvem9saW4uZTEubSA8LSBubHMoCiAgTEwzLmZvcm11bGEsIGRhdGEgPSB2aW5jbG96b2xpbiwgc3Vic2V0ID0gZXhwZXIgPT0gMTA1MDksICAKICBzdGFydCA9IGxpc3QoYiA9IDEsIGQgPSAxMDAwLCBlID0gMC4yNikpCnZpbmNsb3pvbGluLmUyLm0gPC0gbmxzKAogIExMMy5mb3JtdWxhLCBkYXRhID0gdmluY2xvem9saW4sIHN1YnNldCA9IGV4cGVyID09IDEwODIxLCAgCiAgc3RhcnQgPSBsaXN0KGIgPSAxLCBkID0gMTAwMCwgZSA9IDAuMjYpKQp2aW5jbG96b2xpbi5lMy5tIDwtIG5scygKICBMTDMuZm9ybXVsYSwgZGF0YSA9IHZpbmNsb3pvbGluLCBzdWJzZXQgPSBleHBlciA9PSAxMDgyOCwgIAogIHN0YXJ0ID0gbGlzdChiID0gMSwgZCA9IDEwMDAsIGUgPSAwLjI2KSkKdmluY2xvem9saW4uZTQubSA8LSBubHMoCiAgTEwzLmZvcm11bGEsIGRhdGEgPSB2aW5jbG96b2xpbiwgc3Vic2V0ID0gZXhwZXIgPT0gMTA5MDQsICAKICBzdGFydCA9IGxpc3QoYiA9IDEsIGQgPSAyNzAwLCBlID0gMC4wMykpCnZpbmNsb3pvbGluLmU1Lm0gPC0gbmxzKAogIExMMy5mb3JtdWxhLCBkYXRhID0gdmluY2xvem9saW4sIHN1YnNldCA9IGV4cGVyID09IDExMDIzLCAgCiAgc3RhcnQgPSBsaXN0KGIgPSAxLCBkID0gMTAwMCwgZSA9IDAuMjYpKQp2aW5jbG96b2xpbi5lNi5tIDwtIG5scygKICBMTDMuZm9ybXVsYSwgZGF0YSA9IHZpbmNsb3pvbGluLCBzdWJzZXQgPSBleHBlciA9PSAxMTEwNiwgIAogIHN0YXJ0ID0gbGlzdChiID0gMC41LCBkID0gMjYwMCwgZSA9IDAuMDIpKQpgYGAKCgpgYGB7ciwgZmlnLndpZHRoPTMuMywgZmlnLmhlaWdodD0zLCB3YXJuaW5nPUZBTFNFfQojIyBGaWcuIDguNC4gCiNwYXIocHR5ID0gInMiKQpwbG90KGVmZmVjdCB+IGNvbmMsIGRhdGEgPSB2aW5jbG96b2xpbiwKICAgICBwY2ggPSBhcy5udW1lcmljKGV4cGVyKSwgCiAgICAgbG9nID0gIngiLAogICAgIHhsaW0gPSBjKDFlLTA0LCAxMCksIAogICAgIHhsYWIgPSBleHByZXNzaW9uKHBhc3RlKCJDb25jZW50cmF0aW9uKCIsIG11LCAiTSkiKSksIAogICAgIHlsYWIgPSAiTHVtaW5lc2NlbmNlKExVKSIpCmNvbmNWZWMgPC0gZXhwKHNlcShsb2coMWUtMDQpLCBsb2coMTApLCBsZW5ndGgub3V0ID0gNTApKQpsaW5lcyhjb25jVmVjLCBwcmVkaWN0KHZpbmNsb3pvbGluLmUxLm0sIGRhdGEuZnJhbWUoY29uYyA9IGNvbmNWZWMpKSwgbHR5ID0gMikKbGluZXMoY29uY1ZlYywgcHJlZGljdCh2aW5jbG96b2xpbi5lMi5tLCBkYXRhLmZyYW1lKGNvbmMgPSBjb25jVmVjKSksIGx0eSA9IDMpCmxpbmVzKGNvbmNWZWMsIHByZWRpY3QodmluY2xvem9saW4uZTMubSwgZGF0YS5mcmFtZShjb25jID0gY29uY1ZlYykpLCBsdHkgPSA0KQpsaW5lcyhjb25jVmVjLCBwcmVkaWN0KHZpbmNsb3pvbGluLmU0Lm0sIGRhdGEuZnJhbWUoY29uYyA9IGNvbmNWZWMpKSwgbHR5ID0gNSkKbGluZXMoY29uY1ZlYywgcHJlZGljdCh2aW5jbG96b2xpbi5lNS5tLCBkYXRhLmZyYW1lKGNvbmMgPSBjb25jVmVjKSksIGx0eSA9IDYpCmxpbmVzKGNvbmNWZWMsIHByZWRpY3QodmluY2xvem9saW4uZTYubSwgZGF0YS5mcmFtZShjb25jID0gY29uY1ZlYykpLCBsdHkgPSAiMzMxMyIpCmBgYAoKCmBgYHtyfQp2aW5jbG96b2xpbi5tMSA8LSBubG1lKAogIGVmZmVjdCB+IGQvKDEgKyBleHAoYiAqIChsb2coY29uYykgLSBsb2coZSkpKSksIAogIGRhdGEgPSB2aW5jbG96b2xpbiwKICBmaXhlZCA9IGxpc3QoYiB+IDEsIGQgfiAxLCBlIH4gMSksIAogIHJhbmRvbSA9IGQgfiAxIHwgZXhwZXIsIAogIHN0YXJ0ID0gYygxLCAxMDAwLCAwLjEpKQpzdW1tYXJ5KHZpbmNsb3pvbGluLm0xKQpgYGAKCgpgYGB7ciwgZmlnLndpZHRoPTMuMywgZmlnLmhlaWdodD0zLCB3YXJuaW5nPUZBTFNFfQojIyBGaWcuIDguNS4KI3BhcihwdHkgPSAicyIpCnBsb3QoZWZmZWN0IH4gY29uYywgZGF0YSA9IHZpbmNsb3pvbGluLAogICAgIHBjaCA9IGFzLm51bWVyaWMoZXhwZXIpLCAKICAgICBsb2cgPSAieCIsCiAgICAgeGxpbSA9IGMoMWUtMDQsIDEwKSwgCiAgICAgeGxhYiA9IGV4cHJlc3Npb24ocGFzdGUoIkNvbmNlbnRyYXRpb24oIiwgbXUsICJNKSIpKSwgCiAgICAgeWxhYiA9ICJMdW1pbmVzY2VuY2UoTFUpIikKY29uY1ZlYyA8LSBleHAoc2VxKGxvZygxZS0wNCksIGxvZygxMCksIGxlbmd0aC5vdXQgPSA1MCkpCmxpbmVzKGNvbmNWZWMsIHByZWRpY3QodmluY2xvem9saW4uZTEubSwgZGF0YS5mcmFtZShjb25jID0gY29uY1ZlYykpLCBsdHkgPSAyKQpsaW5lcyhjb25jVmVjLCBwcmVkaWN0KHZpbmNsb3pvbGluLmUyLm0sIGRhdGEuZnJhbWUoY29uYyA9IGNvbmNWZWMpKSwgbHR5ID0gMykKbGluZXMoY29uY1ZlYywgcHJlZGljdCh2aW5jbG96b2xpbi5lMy5tLCBkYXRhLmZyYW1lKGNvbmMgPSBjb25jVmVjKSksIGx0eSA9IDQpCmxpbmVzKGNvbmNWZWMsIHByZWRpY3QodmluY2xvem9saW4uZTQubSwgZGF0YS5mcmFtZShjb25jID0gY29uY1ZlYykpLCBsdHkgPSA1KQpsaW5lcyhjb25jVmVjLCBwcmVkaWN0KHZpbmNsb3pvbGluLmU1Lm0sIGRhdGEuZnJhbWUoY29uYyA9IGNvbmNWZWMpKSwgbHR5ID0gNikKbGluZXMoY29uY1ZlYywgcHJlZGljdCh2aW5jbG96b2xpbi5lNi5tLCBkYXRhLmZyYW1lKGNvbmMgPSBjb25jVmVjKSksIGx0eSA9ICIzMzEzIikKbGluZXMoY29uY1ZlYywgCiAgICAgIHByZWRpY3QodmluY2xvem9saW4ubTEsIG5ld2RhdGEgPSBkYXRhLmZyYW1lKGNvbmMgPSBjb25jVmVjKSwgbGV2ZWwgPSAwKSwgCiAgICAgIGx0eSA9IDEsIGx3ZCA9IDMpCmBgYAoKCiMgLS0tLQ==